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Formation of large-scale structures by turbulence 
in rotating planets 

Abstract 

This thesis presents a newly developed theory for the formation and main¬ 
tenance of eddy-driven jets in planetary turbulence. The novelty is that jet 
formation and maintenance is studied as a dynamics of the statistics of the flow 
rather than dynamics of individual realizations. This is pursued using Stochastic 
Structural Stability Theory (S3T) which studies the closed dynamics of the first 
two cumulants of the full statistical state dynamics of the flow by neglecting 
or parameterizing the third and higher-order cumulants. S3T is an analytical, 
predictive and quantitative theory for turbulence that proceeds directly from the 
equations of motion and provides a way of finding turbulent statistical equilibria 
and determining their stability. Instability of the statistics of the flow signifies 
transition of the turbulent regime to a new regime. 

With this statistical closure large-scale structure formation is studied in 
barotropic turbulence on a /3-plane. By studying the dynamics of the statis¬ 
tical state novel phenomena are predicted such as: the instability of homogeneous 
turbulence to jet formation, the establishment of turbulent equilibria, the predic¬ 
tion of multiple turbulent equilibria, jet merging bifurcations, and the existence 
of latent jets. Although these phenomena cannot be predicted by analysis of the 
dynamics of single realizations of the flow, it is demonstrated that the predictions 
of the statistical theory are reflected in individual realizations of the flow. 

It is further demonstrated that at analytically predicted critical parameter 
values the homogeneous turbulent state undergoes a bifurcation and becomes 
inhomogeneous with the emergence of large-scale zonal and/or non-zonal flows. 
The mechanisms by which the turbulent Reynolds stresses organize to reinforce 
infinitesimal mean flow inhomogeneities, thus leading to this statistical state insta¬ 
bility, are extensively studied for various regimes of parameter values (planetary 
vorticity gradient, dissipation rate and turbulent energy injection rate) and it is 
shown that for small and modest values of planetary vorticity gradient, /3, the 
upgradient fluxes responsible for the formation and maintenance of large-scale 
structure are induced by the Orr mechanism, while for large /3 by resonant wave 
triads. It is demonstrated that the S3T instabilities equilibrate to finite amplitude 
jets, in agreement with the jets that develop in individual simulations. The 
relation between the formation of large-scale structure through modulational 
instability and the S3T instability of the homogeneous turbulent state is also 
investigated and it is shown that the modulational instability results are subsumed 
by the S3T results. 

The study of the S3T stability of inhomogeneous turbulent jet equilibria is also 
presented and the relation with the phenomenon of jet merging is investigated. 
Methods for finding inhomogeneous statistical turbulent equilibria and also for 
studying their stability are developed. 
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Introduction 



Figure 1.1: The atmospheres of planets Jupiter and Earth are turbulent. The turbulence is 
anisotropic and inhomogeneous and the kinetic energy of the flow is concentrated in large-scale 
zonal or wavy jets and large-scale vortices which persist in the flow enhancing the large-scale 
long time range predictability of the flow. Credits: NASA/JPL and NASA/GSFC. 


1.1 Jets on Earth and Jupiter 

Turbulent atmospheric flows in rotating planets are observed to self-organize into 
large-scale structures. These structures vary at a time scale much larger compared 
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to the turbulent eddy motions with which they coexist. Prominent characteristic 
examples are the Earth’s subtropical and polar jet streams or the zonal winds 
in Jupiter and its Great Red Spot (see Fig. 1.1). Changes in the structure or 
the position of the Earth’s jet streams may induce dramatic changes in regional 
weather patterns. Recent such examples are the 2003 and 2010 European heat 
waves and the 2013-14 North American cold wave that were caused by shifts in 
the position of the jets. 

Jet streams (or jets) are strong and narrow quasi-rectilinear air currents found in 
the atmospheres of some planets. In the Earth there are two jets in each hemisphere 
that flow eastwards. The typical structure of the mean zonal winds over the Pacific 
ocean reveals the double jet structure in each hemisphere (cf. Figs. 1.2a,b). In a 
frame rotating with the Earth the jets have typical speeds of 40 m s -1 (145 km h -1 ) 
and may reach 70 ms -1 (250 kmh -1 ) in the winter of each hemisphere. The wind 
maximum of the subtropical jet is located at around 30°N/S and at a height of 
10-16 km (or at a pressure of 200 mb). The polar jet is located at 40°-60°N/S and 
at a height of about 10 km (300 mb). The weaker subtropical jet is much more 
axisymmetric while the stronger polar jet has a pronounced slowly translating 
non-zonal wave component, especially in the Northern Hemisphere, as shown in 
Fig. 1.1, with the jet maxima distributed over an annular region as depicted in the 
schematic of Riehl (1962) in Fig. 1.2c. Due to its spatial and temporal variation 
the polar jet stream does not appear as a prominent feature in plots of the annual 
mean zonal velocity. For example, in Fig. 1.3a only the subtropical jets appear. 

While little information is known about the vertical structure of the winds 
on Jupiter 1 there is a lot of information about the latitudinal structure of the 
winds at cloud level (about 700 mb) obtained from involved cloud tracking 
techniques. These measurements revealed the existence of an alternating jet 
structure consisting of 15 eastward and 15 westward jets located at the latitudes 
that separate the colored belts of the planet (cf. Fig. 1.4a,b). Near the equator the 
speed of the eastward zonal jet exceeds by 120 ms -1 (430 kmh -1 ) the rotational 
speed of the planet, indicating that the Jovian atmosphere has an appreciable 
superrotation at the equator. Moreover, the jets of Jupiter seem to vary very 
slowly despite being embedded in strong turbulent flow. This remarkable fact 
was discovered when the wind measurements made by Voyager 2 and the Cassini 

1 Preliminary evidence suggests that the jets increase below the clouds (Atkinson et ah, 1997). 
Definitive answers are expected from analysis of the gravitometric measurements that will be 
collected by space probe Juno (Kaspi, 2013; Read, 2013). 
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Figure 1.2: Earth’s jet streams: Shown are pressure (mb)-latitude sections of the mean zonal 
wind (ms -1 ) and divergent circulation averaged over (a) the west Pacific (120°E-170°E) and (b) 
the east Pacific (130°W-180°W). Shown are time averages for May 2014. Divergent circulation is 
represented by vectors of combined pressure vertical velocity and the divergent component of the 
meridional wind.) Red (blue) shading and solid (dashed) contours denote eastward (westward) 
zonal wind. The eastward wind maxima at around 30°N/S correspond to the subtropical jet 
streams while the maxima at higher latitudes to the eddy-driven polar jet streams. (Credit: 
NWW, Climate Prediction Center.) (c) The mean position of the subtropical jet (thick solid 
line) and the region (shaded) of principal polar jet stream activity for the northern hemisphere 
(Taken from Riehl (1962)). 



Figure 1.3: Zonal mean cross sections of (a) the zonal wind component, u , (ms -1 ), (b) the 
meridional wind component, v, (ms -1 ) and the inferred flow of atmospheric mass (10 10 kgs -1 ) 
for annual-mean and zonally averaged conditions. (Brackets denote time average.) (Taken from 
Peixoto and Oort (1984).) 
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Figure 1.4: Jupiter’s jets: (a) A cylindrical projections of Jupiter from Cassini images 
(Credits: NASA/JPL/Space Science Institute), (b) The zonal wind speeds measured by Limaye 
(1986) on Voyager 2 images taken on 1979 (red curve) and by Porco et al. (2003) on Cassini 
images taken on 2000 (black curve). (Taken from Porco et al. (2003).) (c) A zoom at the 
eastward jet found at 24°N. (Taken from Sanchez-Lavega et al. (2008).) 



Figure 1.5: The Intertropical Convergence Zone (ITCZ) over the eastern Pacific on 12 
July 2000 becomes visible by the nearly zonal cloud band. The major convective activity that 
drives the upwelling of the Hadley cell has been organized to occur in this narrow zone. Credit: 
NASA/GOES. 


space probes were analyzed and were found to produce nearly identical winds, 
although 20 years intervened between the measurements. The respective wind 
measurements made by the two probes are shown in Fig. 1.4b. Also, the jets 
have a very special shape: the jet maxima (the superrotating flow) are very 
pointed while the westward jets (the subrotating flows) are weaker and blunted 
(cf. Fig. 1.4c). This thesis will provide an explanation for both the constancy of 
the zonal winds and their shape. 

In the Earth, the subtropical jet stream is driven by the large-scale merid¬ 
ional circulation (see Fig. 1.3c) initiated by the enhanced convective activity 
which is concentrated in a narrow zonal band near the equator, called the ITCZ 
(Intertropical Convergence Zone) (cf. Fig. 1.5). This axisymmetric circulation 
produces a nearly angular momentum conserving zonal flow by transferring angu- 


4 







lar momentum from the deep tropics to the poleward upper part of the Hadley 
cell, where the subtropical jet maximum is located (Fig. 1.3c). The polar jet is 
maintained from the momentum convergence of the turbulent motions, which 
themselves owe their existence to the very jet they support. It should be clari¬ 
fied that the turbulent motions responsible for the maintenance of the polar jet 
is the midlatitude turbulence, with typical length scales of 1000 km and time 
scales of a week. This large-scale atmospheric turbulence is often referred to 
in the meteorological literature as the macroturbulence (Held, 1999; Schneider 
and Walker, 2006). While the dynamics of the subtropical jet has been fully 
elucidated (Schneider and Lindzen, 1977; Held and Hou, 1980; Lindzen and Hou, 
1988), the theory for the formation and maintenance of the eddy-driven jets is 
far from complete. This thesis will present a new theory for the formation and 
maintenance of eddy-driven jets in planetary turbulence. 

The idea that smaller scale turbulence transfers momentum and maintains larger 
scale flows, thus fluxing momentum upgradiently, is provocative and has been 
called in the literature, equally provocatively, a “negative viscosity” phenomenon 
(Starr, 1953). This idea, expressed in this manner, seems to violate the natural 
entropic tendency of physical systems towards states of greater disorder, which 
is consistent with the usual downgradient action of diffusion: for example, as 
Reynolds (1883) has shown, high density ink is spread and homogenized by 
turbulent flow. Interestingly, as we will discuss shortly, it has been shown that the 
fine-grain maximum entropy states in barotropic flows correspond to macrostates 
with large-scale jets and vortices (Miller, 1990; Robert and Sommeria, 1991; 
Bouchet and Venaille, 2012). 

Jeffreys (1926) was the first to propose that large-scale atmospheric circulation 
is eddy-driven. Until then people were trying to obtain understanding of the 
general circulation of the atmosphere neglecting the non-axisymmetric dynamics 
of the flow as well as the effect of the mean quadratic eddy stresses (the divergence 
of the Reynolds stresses) on the mean axisymmetric circulation. Jeffreys demon¬ 
strated that an exclusively axisymmetric point of view is inadequate by analyzing 
the zonally averaged momentum balance of a whole column of air located at 
the midlatitudes with westerly (or eastward) winds at the ground. Since the 
momentum of the whole column (about 10 5 N s per m 2 ) is lost at the ground at 
the rate of 0.1 Nm -2 (this is the surface drag), it can be shown that it cannot be 
replenished by the flux of momentum by the observed mean axisymmetric circu¬ 
lations, and thus he concluded that the surface westerlies had to be maintained 
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by the horizontal convergence of momentum from the non-axisymmetric motions 
(the eddy motions), i.e., he argued that the horizontal momentum divergence 
of the non-axisymmetric motions must be responsible for maintaining the mean 
momentum of the column and also for the westerlies at the ground. 2 The paper 
of Jeffreys introduced the idea that the eddy motions in the atmosphere (the 
cyclones) should not be viewed as unstable perturbations to the axisymmetric 
mean circulation but rather as an integral component for the maintenance of the 
very axisymmetric circulation that gives rise to them (for a historical discussion 
cf. Lorenz (1967)). Jeffreys also stressed for the first time that the horizontal 
eddy barotropic dynamics are responsible for the maintenance of the large-scale 
structure. This point of view was further advanced and given a theoretical basis 
by Rossby and collaborators (1947a, b). In the meanwhile the detailed upper 
air observations, that became available after the Second World War, gave solid 
observational support to the proposition that the Earth’s polar jet is supported 
by barotropic upgradient momentum fluxes. Such modern measurements of the 
momentum flux compiled by Peixoto and Oort (1984) are shown in Fig. 1.6. In 
these plots the zonally averaged momentum flux uv is decomposed as 3 

uv — uv + u'v' , (1.1) 

into the momentum flux due to the axisymmetric motions, uv, and the mean 
momentum flux due to the eddies (the motions that deviate from axisymmetry), 
u'v'. It can be seen in Fig. 1.6 that the eddy momentum fluxes are larger than the 
momentum fluxes from the meridional circulations and dominate at the latitudes 
of the polar jet. There is a large eddy momentum flux convergence around 55°N/S 
and at 300 mb, precisely at the location of the polar jet, and at this location the 
momentum fluxes from the axisymmetric motions is negligible, showing that the 
polar jet is maintained by eddy momentum flux convergence. The second region 
of convergence occurs at 30°N/S where the subtropical jet is located (cf. Fig. 1.6b). 
The dominant flux convergence in this case is from the axisymmetric motions 
indicating that the subtropical jet is a result of an axisymmetric circulation. 

2 At the time of Jeffreys the existence of the polar jet was not known and in his 1926 paper 
there is no mention of upper-level jets. Evidence of strong upper-level winds was obtained from 
kite measurements first in Japan in the 1920’s, but the existence of the polar jet as a global 
feature of the climatology was established from aircraft measurements during the Second World 
War. The term “jet stream” was coined by Rossby in 1942 (Lewis, 2003). 

3 u is the zonal velocity with eastward flow being u > 0, v is meridional velocity with northward 
flow corresponding to v > 0 and overbar denotes zonal average. 
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Figure 1.6: Zonal mean cross section of zonal mean of northward flux of zonal momentum 
(cf. (1.1)) by (a) the eddies, [u'v'] (m 2 s -2 ) and (b) the mean meridional circulation, [uv\ 
(m 2 s -2 ). Brackets denote time average for a period of a year. It is clear that the dominant 
contribution comes from the eddies. Note the eddy momentum flux convergence occurs in the 
location of the polar jet while mean meridional momentum flux converges at the location of the 
subtropical jet. (Taken from Peixoto and Oort (1984).) 


The Jovian jets are eddy-driven, like the Earth’s polar jets. This was confirmed 
through systematic and repeated analysis of the turbulent velocity fields at cloud 
level (Ingersoll et al., 1981; Ingersoll, 1990; Salyk et al., 2006; Galperin et al., 
2014). Ingersoll and coworkers demonstrated the upgradient action of the eddy 
momentum fluxes by plotting the eddy momentum flux, u'v ', together with the 
shear of the mean flow, du/dy, as a function of latitude. They noted that these 
two quantities are positively correlated and to a great degree of accuracy satisfy 
the linear law, 


du 

'u- ? 
d y 


( 1 . 2 ) 


with k ~ 10 6 m 2 s _1 , as can be seen in Fig. 1.7. That n > 0 implies the remarkable 
fact that the momentum fluxes on Jupiter act anti-diffusively, since the rate of 
change of zonal momentum (disregarding dissipation) obeys 


du 

dt 



d 2 u 


(1.3) 


which is a diffusion equation with negative diffusion coefficient. 


1.2 Theories for jet formation and current understanding 

Since atmospheric motions on Earth are confined in a thin shell 10 000 km in the 
horizontal and 10 km in the vertical (the mean depth of the troposphere, where 
most of the mass of the atmosphere is located) the motions are quasi-horizontal 
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Figure 1.7: The structure of the eddy momentum flux, [u'v'] (dots), and the zonal mean 
flow shear, [du/dy] (solid), as a function of latitude, [u'v'] error bars correspond to 2 standard 
deviations from the mean. There is a distinct positive correlation between the two curves. (Taken 
from Salyk et al. (2006).) 

and one would expect that barotropic dynamics, which involve only the height- 
averaged flow fields, would be sufficient to describe the dynamics of the Earth’s 
atmosphere. However, this is erroneous: the vertical shear of the mean zonal wind 
(i.e. the derivative of the zonal wind with height), which is associated with the 
temperature contrast between the Equator and the Pole, gives rise to powerful 
baroclinic growth processes that produce the cyclones, which are the “atoms” of 
atmospheric turbulence responsible for the transport of heat and momentum in 
the atmosphere. 

Interestingly, the cyclones grow first through the baroclinic process of drawing 
potential energy from the mean flow and transporting in this way heat to the poles, 
and then they assume a nearly barotropic structure (i.e. height independent). 
The barotropic dynamics responsible for the redistribution of momentum in the 
upper troposphere and the formation of jets is also referred to as the “barotropic 
governor”, because this mechanism of barotropic exchange that forms the jets is 
the very mechanism that alleviates the instability of the atmospheric flow and 
maintains the atmosphere in a state of baroclinic neutrality (Ioannou and Lindzen, 
1986; James, 1987; Lindzen, 1993; Roe and Lindzen, 1996). This duality in the 
behavior of the baroclinically growing disturbances simplifies the study of jet 
formation in baroclinic atmospheres. It allows us to consider that the atmospheric 
dynamics fall into two manifolds: the slow barotropic manifold that controls the 
formation and evolution of jets in the upper troposphere, and the faster manifold 
of baroclinic processes that provides the necessary excitation of the barotropic 
manifold to maintain it in a turbulent state. A confirmation of this point of view 
has been given by DelSole (2001), who by considering that the upper troposphere 









is governed by barotropic dynamics excited by the baroclinic activity from lower 
levels demonstrated that the structure of the momentum fluxes responsible for 
maintaining the upper level jets could be accurately captured. As a result, in 
this thesis we will adopt the traditional view and study the formation of jets and 
other large-scale structure both in the Earth and in Jupiter within the context of 
barotropic dynamics. 

This barotropic two-dimensional framework has been adopted by most re¬ 
searchers that investigated jet formation in Jupiter and the outer planets starting 
with Williams (1978) and more recently with Nozawa and Yoden (1997) and 
Huang and Robinson (1998). Other authors investigated the formation of jets on 
Jupiter in the primitive equation extension of the quasi-geostrophic barotropic 
dynamics by modeling the Jovian atmosphere as a shallow-water fluid; but also in 
these studies jet formation proceeded as in the purely two-dimensional barotropic 
models (Cho and Polvani, 1996a,b; Scott and Polvani, 2007). That these dynamics 
can produce jet formation has been also demonstrated experimentally by Read 
et al. (2004) in the Coriolis rotating tank in Grenoble and by Espa et al. (2010) in 
Rome. We adopt the simplest framework and study jet formation in Jupiter and 
the outer planets in the context of barotropic dynamics that are maintained in a 
turbulent state by external excitation. The excitation represents the introduction 
of vorticity at cloud level from convective motions induced by the heating source 
in the interior planet. The typical scale of vorticity injection is 1000 km (Little 
et ah, 1999; Gierasch et al., 2000). 

We now review the main theories that have been advanced for understanding 
the formation of jets in turbulence. These theories can be distinguished as those 
that arise from turbulence theory and are generally phenomenological, and those 
that consider that the flow perturbations are coherent, like a wave, and study the 
interaction of this coherent eddy field with the mean flow. The latter theories 
will be referred to as wave-mean flow interaction theories and are generally more 
deductive. 

1.2.1 Turbulence theories 

Jet formation in turbulence theory is viewed as a consequence of a cascade of 
energy from small scales to large scales. This type of cascade is called “inverse” 
and is the opposite of direct cascades that characteristically operate at small 
scales in homogeneous isotropic 3D turbulence transferring energy from large 
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scales to small scales where it is dissipated. That turbulence confined on a 
plane, like the barotropic turbulence that we will study, supports an inverse 
cascade in energy was first proposed by Fj0rtoft (1953) who argued that this was 
consequence of the two dimensionality of the flow which implies in the inviscid 
limit the conservation of the total kinetic energy of the flow, E — f ||u| 2 dA, 
as in 3D, but also the conservation of the vorticity of every particle in the flow, 
which leads to an infinite set of integral invariants. As a result, on the x-y plane 
the material conservation of the vorticity ( = (V x u) z = d x v — d y u , implies that 
all integrals over the whole area of the fluid of the form f E(() dA, with F any 
integrable function, are conserved. Enstrophy, defined as Z = / ^£ 2 d A, is the 
invariant that is usually considered out of this hierarchy of conserved quantities 
and Fjprtoft considered the constraint imposed on the spectral evolution of the 
flow by the simultaneous conservation of energy and enstrophy. Expressing the 2D 
incompressible flow field through a streamfunction ^asu= (u, v ) = (—dy'ip, d x ^), 
implies that ( — (d 2 x + <9 2 y )^, and expanding the streamfunction in Fourier as, 
if; = (2?r)~ 2 f d 2 k c lk ' x , we have that energy and enstrophy are respectively 
given as E — (2tt)~ 3 f d 2 k |k| 2 1| 2 and Z — (2ir)~ 3 f d 2 k |k| 4 1| 2 . This means 
that the energy and enstrophy spectral densities that correspond to wavenumber 
k, £(k) and Z( k) respectively, are related through Z{ k) = k 2 £( k), where k — |k|. 
Fjprtoft stated that conservation of energy and enstrophy in 2D flows constrains 
the energy exchanges between scales in such a manner so that if enstrophy moves 
to smaller scales energy must move to larger scales. 

However, these energy exchanges in the unforced, inviscid limit are reversible 
in time and as a result no systematic direction of energy or enstrophy flow can 
be deduced from these arguments (Salmon, 1998; Tung and Orlando, 2003b). In 
irreversible forced-dissipative systems, it can be argued that energy should move 
to large scales and enstrophy to small scales, as Fjprtoft envisioned, namely from 
the scale the energy or enstrophy is being injected towards the scale that each is 
dissipated. Kraichnan (1967) provided such a refinement of Fjprtoft’s argument, 
which was further refined by Eyink (1996), by showing that energy and enstrophy 
conservation imply an inverse energy cascade since if energy and enstrophy are 
injected at a scale kf energy must be dissipated at a larger scale k r < kf and 
enstrophy at a smaller scale, k v > kf. A The physical mechanism that decreases 

4 Kraichnan’s argument: Assume that energy and enstrophy are injected at a scale kf at 
rates Sf and 77 / = k 2 f Sf, and that they are being dissipated at two distinct scales: a larger scale 
k r (with k r < kf) at rates e r and 7 y r , and a smaller scale k v > kf at rates £ v and r/ v and that 
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Figure 1.8: Evolution of the vorticity field in forced-dissipative flow on a doubly periodic 
domain, 2i r x 2i r. Energy is being injected in the flow at rate Sf and at length scale 2ix/kf with 
kf = 40, and dissipation is done at small scales with 8th-order hyperdiffusion with coefficient 
v — 7 x 10 -35 . (a) Typical structure of the forcing field, (b) The vorticity field at time t/r = 1265, 
with r = (e /&/) -1//3 • Vorticity patches tend to merge creating large-scale vortices which contain 
most of the energy in the flow. As vortices are advected by the flow characteristic vortex filaments 
are created transferring enstrophy to smaller scales, (c) Velocity held at time t/r — 1265. For 
the simulation a pseudospectral code was used at a 256 X 256 resolution with an exponential 
filter acting on wavenumbers k > |/c ma x, where k max is the maximum resolved scale. 

the mean scale of the enstrophy of the flow is the stretching of the vorticity as the 
flow evolves. This is shown in a simulation of forced-dissipative 2D turbulence in 
Fig. 1.8. Vortex filaments form increasing the vorticity gradients and transferring 
enstrophy to smaller scales and, as has been argued by Fjprtoft, energy to larger 
scales in the form of large vortices. 

Kraichnan (1967); Leith (1968); Batchelor (1969) (ofter abbreviated as KBL) 
suggested that conservation of energy and enstrophy in 2D turbulence results 
in the formation of two distinct inertial ranges: a range in which energy is 
transferred upwards to larger scales and a range in which enstrophy is transferred 
to smaller scales. Using Kolmogorov type non-dimensional arguments, which 
assume that these ranges are homogeneous, isotropic and self-similar, they showed 
that the energy density in the energy transferring inertial range should follow 
the power law £(k) = where e r is the energy transfer rate and C 

there is almost no dissipation for k r < k < k v . At a statistical steady state we expect from 
conservation of energy and enstrophy that, £f = e r + e v and rjf = rj r + rj v , from which we obtain 

£ y _ 1 — ( k r /kf ) r h_ — ( 1 ~ (k r /kf) /-< a\ 

Er~ {k v /k f ) 2 - 1 ’ rjr~\kf) Ur/ (k v /kf) 2 ~ l' 

Note that because all £ and r] are positive such a steady state is possible only if kf satisfies 
k r < kf < k v , as it was assumed from the start . Moreover, if k v kf then £ v /e r <C 1, which 
means that that most of energy flux occurs at large scales k r . Also, if k r <C kf then f] v /r]r 1 
implying that there is no enstrophy flux at large scale and all the enstrophy moves to small scale 
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a dimensionless universal constant, while the energy density spectrum in the 
enstrophy inertial subrange follows the power law £(k) — C'r]v^k~ 3 , where r] v is 
the enstrophy transfer rate and C' a different dimensionless constant. Evidence 
of this scalings has been verified in numerical simulations, as shown in Fig. 1.9a. 
Experiments with flowing soap films in which turbulence is excited by grids (arrays 
of cylinders) lining the walls of the flow channel, provide a physical occurrence of 
forced-dissipative 2D turbulence and confirm the predictions of KBL for the two 
inertial ranges, as seen in Fig. 1.9b. 

All these arguments however, are based on homogeneous and isotropic 2D 
turbulence. It may be the case that large-scale processes in the atmosphere, 
including jet formation, are essentially two dimensional, but overall the atmosphere 
is neither isotropic nor homogeneous. Anisotropy in the atmosphere is induced by 
the Earth’s rotation which distinguishes the zonal from the meridional direction 
and also by the temperature difference between the Equator and the Pole, while 
homogeneity is broken by the large-scale jets. Therefore one should use the 
classical KBL 2D arguments with caution when trying to explain atmospheric 
motions. Also, the observed atmospheric spectrum of the winds in the upper 
troposphere, contrary to what the classical 2D picture would expect, shows 
the fc -5 / 3 dependance on the short-wave side of the spectrum at scales ranging 
from 600 km down to 2 km, at the so called “mesoscales”, as seen in Fig. 1.9c. 
Coincidentally, the direct energy cascade that is typically found at homogeneous 
isotropic 3D turbulence and is responsible for transferring energy to the smaller 
scales where is dissipated, also presents a fc -5 / 3 dependance. However, the vertical 
extend of the troposphere limits 3D effects in the atmosphere to appear only 
at most at scales of 1-2 km and therefore the classical dimensional arguments 
that are offered to account for the fc -5 / 3 inertial range scaling cannot apply 
to the mesoscales. A concrete explanation for the atmospheric spectrum is an 
open and challenging problem, which will not be addressed in this thesis. It is 
interesting to note that if the atmosphere is represented crudely as a two-layer 
fluid, which is one step of an approximation higher than the one adopted in this 
thesis, the atmospheric spectrum of Nastrom and Gage (1985) is obtained (Tung 
and Orlando, 2003a). 

There is an additional problem with the predictions of the KBL theory in 
planetary turbulence. While the existence of the inverse energy cascade could 
provide an explanation for the emergence of the observed large-scale flows in 
planetary turbulence, the theory predicts that the inverse cascade will lead to the 
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Figure 1.9: (a) Spectral energy density E(k) vs total wavenumber obtained by numerical 

simulation of forced-dissipative isotropic two dimensional turbulence, for different resolutions: 
(A) 4096 x 4096 and (B) 32 768 x 32 768. Also shown are the k ~ 5//3 (dashed) and k~ 3 (dotted) 
slopes. (After Boffetta and Musacchio (2010).) (b) Experimental spectral energy density at 
steady state of flowing soap films in which turbulence was excited by grids (arrays of cylinders) 
lining the walls of the flow channel, showing clear evidence of the dual inverse energy and 
forward enstrophy cascade. The k ~ 5//3 and k~ 3 slopes are also shown. (After Rutgers (1998).) 
(c) Variance power spectra of winds near the tropopause from commercial aircraft data (by 
NASA GASP). The spectrum for meridional wind is shifted one decade to the right so it is 
visible. (After Nastrom and Gage (1985).) 
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formation of a large-scale condensate, as large as the geometry allows, as shown 
in Figs. 1.8b,c. The structures that emerge in planetary turbulence however 
are usually not at the largest scale of the flow and moreover they have a very 
particular structure (see for example Fig. 1.4c). As a result if we are to provide a 
theory for the emergence and maintenance of the large-scale structure in planetary 
turbulence, we have to go beyond the classical KBL theory and consider the 
implications of anisotropy and even inhomogeneity in 2D turbulence. 

1.2.2 Rossby waves 

Flows at rotating planetary atmospheres are anisotropic due to the preferential 
direction imposed by the rotation. This has important implications to planetary 
motions because it leads to a new class of exact solutions of the equations of 
motions that were discovered by Rossby and are referred to as Rossby waves 
(Rossby, 1939). Consider a fluid on the nearly spherical surface of a rotating 
planet at angular velocity fi, where the magnitude fl = |fi| is the rate of rotation 
of the planet (for the Earth D = 7.27 x 10 -5 rads -1 = 2n (24h) -1 , for Jupiter 
ft = 1.76 x 10 -4 rads -1 ~ 2n (10h) -1 ). The velocities of the fluid as measured in 
an inertial frame of reference and as measured in a frame of reference co-rotating 
with the planet are related by u/ = ur + fi x r, where subscripts / and R denote 
quantities measured in the inertial and rotating frame respectively. That the 
velocity of the flow lies predominantly on the plane tangential to surface of the 
sphere implies that the vorticity in the inertial frame, V x u/ ~ z, is normal 
to the surface of the sphere (z being the direction normal to the surface of the 
planet as shown in Fig. 1.10). The magnitude of the vorticity, £/, is equal to 
£/ = (r + /, where (r = (V x u r) • z is the relative vorticity of the fluid and 
/ = (2ft) • z, is the planetary vorticity, which is also referred to as the Coriolis 
parameter being the coefficient that multiplies the velocity in the Coriolis force 
in the equations of motion in a rotating frame. Material conservation of vorticity 
on this planetary surface implies that is conserved, i.e., 

^Dt = D t^ R + ^ = ip 1 + (Cr + /) = 0 5 (l- 5 ) 

where D/D t = dt + u • V is the material derivative that determines the time 
rate of change of flow properties as they move with the fluid, which is frame 
independent when it acts on scalars. From hereafter we will drop subscripts R for 
simplicity and, as it is usually done, all velocities will be considered to be relative 
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Figure 1.10: A /3-plane at latitude Oo approximates a zonal belt of the spherical surface 
centered at Oo with a plane. The coordinate x in that plane corresponds to the zonal direction 
on the sphere (circles of constant latitude 6) and the direction y is taken to correspond to the 
meridional direction on the sphere (circles of constant longitude A). 

to the rotating frame of reference. 

The spherical shape of the planet implies that / = 2 sin 6 and consequently 
u • V/ = /3u, with v the poleward velocity, /3 = (212 cos 6) /a, a the radius of the 
planet and 6 the planet’s latitude (cf. Fig. 1.10). This /3-term in the equations of 
motions introduces a new class of anisotropic exact nonlinear solutions, discovered 
by Rossby (1939), that will be of principal importance in this thesis. Rossby 
further introduced the /3-plane approximation that greatly simplifies the equations 
of motions. Instead of solving for the barotropic dynamics (1.5) on the surface of a 
sphere we approximate the domain as a planar surface tangent to the surface of the 
sphere at latitude 6q rotating at the constant rate flsin^o- The anisotropy due to 
the sphericity is then simply introduced by keeping only the /3-term in the equations 
of motion with /3 = (2Q cos 9q)/ a (at latitude 45°, /3 = 1.61 x 10 -11 m -1 s -1 on 
the Earth and /3 = 3.56 x 10 -12 m -1 s -1 on Jupiter). Following the seminal work 
of Rossby and the multitude of theoretical work that followed that employ the 
/3-plane approximation in the study of midlatitude planetary dynamics, we will 
also adopt in this thesis the /3-plane approximation for studying the formation of 
jets and other large-scale structure in planetary turbulence. 

We present briefly the basic features of the Rossby wave solutions which will 
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be a central to this thesis. Expressing the two dimensional velocity in terms of a 
streamfnnction as u = z x Vf, the barotropic vorticity equation (1.5) takes the 
form 

d t A^ + J (^, A= 0 , (1.6) 

where J(g,h) = ( d x g)(d y h ) — ( d y g)(d x h ) is the Jacobian of functions g and h 
and A g = (d% x + d yy )g the horizontal Laplacian. Equation (1.6) reveals that the 
/3-term supports the wave solutions ^ with wavevector k = ( k x , k y ) 

and the frequency satisfying the anisotropic dispersion relation 


w(k) = 



(1.7) 


which is described geometrically in the manner of Longuet-Higgins (1964) in 
Fig. 1.11. Being anisotropic, wavevectors k that correspond to frequency uj do 
not lie on a circle centered at the origin. Figure 1.11 demonstrates the important 
property that Rossby wave packets have y -phase velocities opposite to their 
y-group velocities. 5 

Remarkably, because these monochromatic Rossby waves satisfy J(^, A-0) = 0, 
they are also nonlinear solutions of (1.6). Moreover, stationary Rossby waves 
with k x — 0 correspond to sinusoidal zonal jets with streamfunction = Ae lkyV 
or more generally any zonal flow with ^ = / A{k y )e lkyV d k y is also nonlinear 
solution of (1.6) (other nonlinearly valid Rossby wave solutions are presented in 
Appendix G). We will demonstrate in this thesis that the emergence of large-scale 
features in turbulent flow in the form of jets and other Rossby waves (or zonons) 
can be traced to the property that exact nonlinear solutions can serve as good 
repositories into which the eddy energy may “condensate”. 


5 For example in Fig. 1.19a, which will be discussed later, the phase lines in the region y > 3.6 
are such so that the group velocity is directed towards the center of the channel. Also the phase 
lines in the region y < 2.6 are configured so that group velocity is also directed towards the 
center. 
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Figure 1.11: (a),(c) Lines of constant phase of a Rossby wave with wavevector k subtending an 
angle 0 — arctan(|/c y |/|fcc|) with the horizontal. (b),(d) The locus of wavectors k that correspond 
to a fixed value of uj > 0 is the circle shown with center at C and radius /3/(2uj). Expressing (1.7) 
as uj — (/3/k) cos 6 we find that the phase velocity is c p = (/ 3/k ) (l,tan0), parallel to k, while 
the group velocity is c G = Vk uj = —(/3/k 2 ) (cos(2$), sin(20)), in the direction W&. The y- 
component of the group velocity of the Rossby waves is always directed in the opposite direction 
from the y-component of the phase velocity (phase and group velocity directions are drawn in 
(a) and (c)). 


1.2.3 Anisotropic turbulence on a beta-plane 

The presence of Rossby waves affects the structure of turbulence in the degree that 
the /3-term dominates the advection term, J(Vb A?/)), in the vorticity equation (1.6). 
The ratio of these terms is \/3 d x ip\/\J (ip, A ip) | = WR oss b y /w t urb = O (/?/( k 2 U ims )), 
where cjRossby = 0(/3/k) is the Rossby wave frequency and c^turb = kU rms is the 
inverse of the shear time associated with nonlinear advection, with U Trns the 
root-mean-square velocity of the flow at scale. It is reasonable to expect that 
when cjRossby c^turb the /3-term barely affects the turbulent dynamics. Since 
cjRossby increases while ur b decreases as k decreases, we expect only the large 
scales to be influenced by /3. The wavenumber k Rh = \J(3/U YYns at which the shear 
time scale is equal to the Rossby wave period separates, according to Rhines, the 
wavenumber space in two regions: a region of wave turbulence in which coherent 
Rossby wave motions are manifest in the flow and nonlinearly interact as waves 
and a region in which wave motions are not discernible, the flow is considered 
incoherent and the nonlinear interactions are no longer constrained to be among 
waves. The scale that separates these two regions is being referred to as the 
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Figure 1.12: The interior of the dumbbell or lazy eight (shaded) is the region in which 
Rossby waves are prevalent. According to Rhines (1975) the excitation of the Rossby waves 
halts and anisotropizes the cascade, energy proceeds as shown towards the k x — 0 axis. 

“Rhines’s scale” and the locus of the wavenumbers that satisfy the requirement 
^Rossby — ^turb is the popular and iconic dumbbell shape of Vallis and Maltrud 
(1993), shown in Fig. 1.12, and which will be encountered in this thesis under a 
different exegesis. Rhines (1975) argued that for scales larger than the Rhines’s 
scale (k < k Rh) the selectivity imposed to the nonlinear interactions by wave 
turbulence, which allows interactions only among waves, retards the inverse energy 
cascade and anisotropizes it. Rhines in this way explained that in anisotropic 
/3-plane turbulence the inverse cascade should not be expected to proceed to the 
largest scale available and is halted at k rr This led to the general prediction 
that the large-scale structure in /3-plane turbulence should have a characteristic 
length scale of the order of l//cRh- 

Still unanswered remains how zonal jets finally emerge. Broadly speaking three 
different approaches attempt to answer this question within Rhines’s phenomenol¬ 
ogy. Rhines (1975); Vallis and Maltrud (1993); Chekhlov et al. (1996); Smith 
and Waleffe (1999) argue that the cascade proceeds through local interactions 
up to the dumbbell where the cascade process is halted and the upscale flow of 
energy is directed to move tangentially along the dumbbell (in the direction of 
the arrows in Fig. 1.12) towards the bottleneck at k x — 0, thereby forming jets. 
McIntyre, Dritschel, Scott and collaborators (Baldwin et al., 2007; Dritschel and 
McIntyre, 2008; Dritschel and Scott, 2011; Scott and Dritschel, 2012) argue that 
jet formation is the inevitable and universal result of irreversible mixing of the 
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potential vorticity (PV) of the flow, q = £ + /o + /3y, that occurs in turbulent 
flows, which tends to wipe out the large-scale PV gradients producing a flow 
that approximately satisfies 'Vq = 0 or the large-scale flow satisfies d xx v = d xy u 
and dy y u — d xy v = /3. They further postulate that because the zonal direction, 
x, is a homogeneous direction the resulting well mixed large-scale flows must 
be independent of x and consequently irreversible mixing in the presence of /3 
produces only mean zonal jets u(y) at large scale with a parabolic profile satisfying 
d 2 u/dy 2 = /3 (as for example the central jet of Fig. 1.13b). Their argument has 
however a further twist: there is an important feedback between Rossby waves 
and turbulence that occurs in the boundary of the dumbbell. When waves are 
strongly present the turbulent mixing is inhibited, while when the turbulence 
is strongly nonlinear turbulent interactions iron out the PV gradient. In the 
presence of a zonal jet, the effective Rossby restoring force (“Rossby elasticity”) 
is not /3 but rather, 


R _ d g _ a ^ 

Peff i P j o ’ 

dy dy A 


(1.8) 


which means that in regions of eastward jet maxima Rossby wave excitation is 
reinforced since /3 e ff > /3 and c^Rossby increases inhibiting mixing, while at westward 
jet maxima the excitation of Rossby waves is not comparatively favored, since 
/3 e ff < /3, and this leads to increased mixing of PV that reduces further the PV 
gradient. These two effects result in inhomogeneous PV mixing in the anisotropic 
y direction producing a staircase PV profile, as in Fig. 1.13a, with regions in 
which the PV gradients have been severely reduced and the flow is retrograde 
and parabolic and regions in which the potential vorticity gradient is very large 
with very sharp prograde jets, as shown in Fig. 1.13b. It is remarked by McIntyre 
that the same interaction mechanism was proposed by Phillips (1972, 1977) in 
order to account for the widespread occurrence of a succession of layers of uniform 
stratification in the stratified ocean. 

The observation that the maintained jets in planetary /3-plane turbulence are 
prograde jets joined with parabolic wind profiles is astute. It gives the shape of 
the 24°N jet on Jupiter, shown in Fig. 1.4c, and of the equatorial retrograde jets 
on Neptune and Uranus (see Fig. 1.14a,b). In numerical steady state simulations 
of an almost inviscid turbulent flow on a doubly periodic /3-plane channel Danilov 
and Gurarie (2004) were able to produce the mean flow shown in Fig. 1.14c, which 
almost exactly conforms to the above specification. With this successful prediction 
it is very tempting to cease further effort and accept that the phenomenon of 
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Figure 1.13: Jet formation through homogenization of PV, # = £ + /. Shown are (a) zonal 
mean PV profiles, q = fo + /3y — du/dy , and (b) the corresponding zonal mean flows, u. A linear 
PV distribution (dashed) is associated with no zonal mean flow. Westward jets are associated 
with regions of homogenized PV. These jets become stronger as the homogenization increases 
(compare the jets associated with the dashed-dot PV partially homogenized distribution with the 
staircase profile). Note that PV mixing theory implies that PV gradient should be everywhere 
non-negative. 





Figure 1.14: (a),(b): Observed zonally averaged zonal winds on (a) Uranus and (b) Neptune at 
cloud level. Data taken from Voyager 2 (circles) and from Hubble Space Telescope measurements 
(squares). Solid lines are empirical fits on the data. The equatorial retrograde jets in both planets 
(up to approximately 25°N/S for Uranus and 35°N/S for Neptune) can be approximated very well 
by a parabola. (Taken from Kaspi et al. (2013).) (c) Numerical simulation of forced-dissipative 
turbulence on a /Uplane with large-scale dissipation in the form of second order hypofriction, 
—kA~ 2 . Shown are the zonal mean vorticity gradient, d£/dy, (thin line) averaged over 100 
time units after the simulation has reached steady state and also zonal velocity, u , (thick line) 
multiplied with 100 to fit the same axes. For the simulation a pseudospectal code at 512 x 512 
resolution is used and turbulence is maintained against dissipation by energy injection in the form 
of isotropic excitation at wavenumber kf. The parameters of the simulation are: fikffn = 49 
and £f kf/K 3 = 7.7 x 10 4 . (Taken from Danilov and Gurarie (2004).) 
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Figure 1.15: (a) A snapshot of the vorticity held of forced-dissipative turbulence on a /3-plane 
with the associated instantaneous zonal mean zonal how, u , also indicated. The zonal mean how 
contains 45% of the total energy of the how. (b) The energy spectrum of the zonal how, Ez 
(red), and the spectrum of the remaining energy, Er (blue). The turbulence is forced at the 
scale kf and Er develops the universal k~ 5//3 spectrum while the large scales develop a sharper 
k~ 4 -k~ 5 spectrum as expected from the near discontinuity of the prograde jets shown in (a). A 
pseudospectal code was used at 512 x 512 resolution with /3/(kfr) m 42 and £f k'j/r 3 = 2.3 x 10 7 , 
with r the coefficient of linear damping. Also plotted are the k~ 4 , k~ 5 and k~ 5//3 slopes (dashed). 
(Taken from Bakas and Ioannou (2015).) 


jet formation has in essence been resolved by the above inhomogeneous mixing 
arguments. However, the above arguments are phenomenological, qualitative, 
descriptive and they do not comprise a deductive theory that proceeds from the 
equations of motion. In this thesis we will present a predictive and quantitative 
theory that proceeds directly from the equations of motion that can account 
for the observations. Moreover, this theory leads to different conclusions about 
the role of the physical processes that lead to jet formation. For example, it 
will be shown that the tendency for jet emergence occurs even in the absence 
of /? (in fact /? may actually retard the tendency for the emergence of jets) and 
that /? is required in order to obtain steady state and hydrodynamically stable 
equilibrated flows, which requires that the PV gradient, dq/dy = (3 — d 2 u/dy 2 , be 
of one sign (here positive) and as result in the retrograde parts of the flow, where 
d 2 u/dy 2 > 0, the equilibrated flow must become parabolic satisfying /? = d 2 u/dy 2 , 
while the prograde jets can become infinitely sharp with no constraint other than 
mean diffusive dissipation. The tendency towards discontinuity of the derivative 
of the mean flow implies that the turbulent spectra at large scales have a fc -4 
power law behavior, which is not far from the observed spectrum, shown 

in Fig. 1.15. Note also that the demand that the mean flow does not violate the 
Rayleigh-Kuo criterion for barotropic instability also sets the scale of the jets to 
be nothing else but the Rhines’s scale (/cr h = ^//3/U) given that the maximum 
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of the mean flow velocity [/, rather than r.m.s. turbulent velocity, must at most 
satisfy the condition k 2 U — /3. More basically, this scale should be expected to 
emerge irrespectively of mechanism because it is the only length scale that can 
be formed from the mean flow U (units LT -1 ) and /3 (units L -1 T _1 ). 

For scales within the dumbbell the /3-term dominates over the advection term 
and the flow is well approximated by a sea of weakly interacting Rossby waves. By 
transforming (1.6) in Fourier space it can be shown that the nonlinear Jacobian 
is transformed into a convolution with the property that Fourier component p 
of the flow interacts with Fourier component q to produce Fourier component k 
only if the wavevectors form a triangle, i.e., satisfy p + q = k. In the wave regime 
however, the interacting Fourier components must produce a wave motion and 
significant response is obtained only when the wavenumbers satisfy additionally 
the resonant condition 6 : 

w(k) = w(p) +w(q) . (1.9) 

The turbulence that results from these resonant interactions among waves is 
referred to as “weakly nonlinear turbulence” or “wave turbulence” (WT) (Za¬ 
kharov, 1965; Zakharov et al., 1992; Hasselmann, 1966, 1967). Balk discovered 
that Rossby waves in the WT regime do not only conserve energy and enstrophy 
in the inviscid limit, but also new independent invariant. Balk with Zakharov and 
Nazarenko have demonstrated that this additional invariant, which they named 
“zonostrophy”, is responsible for the anisotropization of the cascades and leads to 
the emergence of zonal jets (Balk, 1991; Balk et al., 1991; Nazarenko and Quinn, 
2009; Balk and Yoshikawa, 2009). 

Finally, a theory of very different character has been advanced for the emergence 
of zonal jets which is based on the property that basic flows consisting of infinitely 
coherent monochromatic Rossby waves are hydrodynamically unstable to zonal 
jets (Lorenz, 1972; Gill, 1974). This instability is called a “modulational instability” 
(MI) because of its similarity with the Benjamin-Feir instability of surface gravity 
waves (Benjamin, 1967; Benjamin and Feier, 1967; Yuen and Lake, 1980; Zakharov 
and Ostrovsky, 2009) and has recently resurfaced in relation to zonal jet formation 

6 Consider a streamfunction ^ T^q, which is a sum of two monochromatic Rossby waves 
i/j p = and ^q = The advection term in (1.6) is then given as: 

Aip) = -A 2 (p 2 - q 2 )(p x q) • ^ e i{(p+ q )- x -Kp)+udq)]t} • 

In the special case for which cj( p) + a;(q) = o;(p + q) the r.h.s. is proportional to a third Rossby 
wave ^p+q which can then resonantly grow. 
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in planetary turbulence and also in drift-wave turbulence in plasmas, both of 
which are governed by the Charney-Hasegawa-Mima equation, which is formally 
equivalent to the barotropic vorticity equation with finite radius of deformation 
(Connaughton et al., 2010). The MI theory for the emergence of jets departs 
considerably from the cascade theories of jet formation. It does not require 
that jets emerge through a sequence of local interactions in wavenumber space 
transferring energy upscale, but rather jet emerge from the instability of the 
primary Rossby wave to a zonal jet perturbation, an interaction that involves a 
non-local interaction in wavenumber space, i.e., the primary Rossby wave with 
wavenumber k = (k x , k y ) gives energy to the spectrally removed unstable zonal 
jet p = (0,p y ), which is nothing else but a zero frequency Rossby wave. In this 
thesis we will investigate the relation of the MI theory for the emergence of jets 
with the statistical theory that will be studied in this thesis. We will show that 
MI is a special case of the more general instability that occurs in the theory we 
will present (cf. chapter 4). 

1.2.4 Statistical approaches for large-scale structure formation 

Turbulent flows involve enormous complexity and a large number of degrees 
of freedom so it is tempting to describe turbulent flows by statistical methods 
reducing in this way its complexity, in a similar way thermodynamics dramatically 
reduce the complexity of the gas molecules movements in a box while still fully 
describe the gas macrostate. However, turbulent flows are usually far from 
equilibrium and the application of equilibrium statistical mechanics might not 
be possible. Interestingly, due to the enstrophy conservation in 2D flows an 
equilibrium statistical mechanical description of the flow in the inviscid limit is 
feasible. 7 The statistical mechanics of a set of inviscid point vortices in 2D flow 
goes back to Onsager (1949) (cf. Eyink and Sreenivasan (2006)) and the first 
statistical mechanical formulation of 2D unforced inviscid flows was developed by 
Miller (1990) and Robert and Sommeria (1991) (the MSR theory). 8 Their theory 
postulates that the emergent equilibrium structures will be the ones that maximize 
entropy while conserving energy, enstrophy and all the hierarchy of invariants in 
2D. Bouchet and Sommeria (2002) extended the MSR theory to quasi-geostrophic 

7 In 3D flows, there can be energy dissipation due to vortex stretching even in the inviscid 
limit, a phenomenon called “anomalous dissipation” (Onsager, 1949; Kaneda et ah, 2003), not 
allowing 3D flows to reach equilibrium state. 

8 For a historic review of statistical mechanical methods in turbulence refer to Bouchet and 
Venaille (2015). 
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barotropic turbulence and showed that the most probable structures are zonal 
jets or large-scale vortices (for a review see Bouchet and Venaille (2012)). 

However, planetary flows are both strongly forced and dissipated and therefore 
out of equilibrium. A non-equilibrium statistical approach is more suitable 
for the description of the statistical dynamics of the turbulent state. Such a 
non-equilibrium statistical mechanical theory has been advanced by Farrell and 
Ioannou (2003) and will be discussed in this thesis. 

1.2.5 Wave-mean flow interaction theories 

It has been known for a long time that waves in a material medium can interact 
with the medium and form mean flows. The acoustic streaming experiments of 
Rayleigh demonstrate this phenomenon (cf. Rayleigh (1896); Lighthill (1978)). 
In acoustic streaming strong jet-like winds are generated by powerful ultrasound 
sources. Acoustic streaming results from the divergence of Reynolds stresses 
induced by the acoustic waves as they dissipate. As with acoustic streaming, jets 
in planetary atmospheres can emerge from Rossby wave streaming in the presence 
of dissipation. This is the basis of the wave-mean flow interaction theories for 
the emergence of jets in the atmospheres. For this mechanism to work the region 
of excitation of the waves and the region of dissipation of the waves should be 
separated. In this case prograde flows emerge in the excitation region while 
retrograde flows emerge in the regions of dissipation. 9 In the Earth the source 
of the equivalent barotropic planetary waves in the upper troposphere, following 
Kuo (1951); Hoskins (1983); Held and Hoskins (1985) and as discussed earlier, 
is the baroclinic activity in the lower troposphere. The equivalent barotropic 
Rossby waves are radiated way to the North and South of the source region where 
they eventually dissipate maintaining the upper-level polar jets (for a model of 
this see DelSole (2001)). The presence of dissipation the emergence of large-scale 
mean flows at steady state since in this case the wave-mean flow non-interaction 
theorem (Eliassen and Palm, 1961; Charney and Drazin, 1961; Andrews and 
McIntyre, 1976; Boyd, 1976b,a) does not hold. (The dissipation of the flow is 
mainly due to Ekman spin-down and additionally to breaking of the waves at 
the “equatorial surf zone” at the critical layers in the Equator-wards flank of the 
midlatitude jet.) 

9 Because of momentum conservation the integrated mean flow acceleration induced by the 
waves integrates to zero. 
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zonal velocity 


Figure 1.16: Schematic explaining the emergence of the upper-level eddy-driven tropospheric 
jet in the Earth’s atmosphere according to Kuo (1951); Held and Hoskins (1985). Equivalently 
barotropic waves are excited by baroclinic processes at the baroclinically active latitudes and as 
they propagate way from the region converge prograde momentum into the region giving rise to 
westerly flows. The waves are dissipated far from the region of excitation forming easterly flows 
at these latitudes. 


To demonstrate this mechanism consider a Rossby wave source region. Rossby 
waves that radiate away from this region converge wave momentum into the 
excitation region, i.e., d y {u'v') < 0, inducing a positive mean flow acceleration 
in this region. This is because Rossby waves with positive group velocities 
propagating to the North (cf. Fig. 1.11a) have u'v' < 0, and Rossby waves with 
negative group velocities propagating to the South (cf. Fig. 1.11c) have positive 
u'v' >0. As a result d y (u'v') < 0 in the stirring region and because the zonal 
mean flow is governed by dtu — —d y (u'v'), a positive mean flow acceleration 
occurs in the stirring region. In the regions of dissipation momentum divergence 
leads to the emergence of negative flows, as shown in Fig. 1.16. If the wave 
excitation is statistically stationary the momentum convergence, —d y (u'v'), will 
also be statistically stationary and the eastward mean flow will grow in the 
mean linearly at a rate proportional to the energy input power, Sf, as shown in 
Fig. 1.17. 10 In conclusion: this wave-mean flow mechanism for the emergence of 
flows predicts linear mean growth of the jet in the regions of stirring and requires 
a localized forcing region and a propagation mechanism (here guaranteed because 
of the positive PV gradient due to the existence of /3) in order for the waves to 
dissipate away from the source region. As a corollary: if the forcing is distributed 
homogeneously and the dissipation coefficients are constant (with no preferred 
dissipation regions) then this mechanism cannot induce any mean flows. 

10 This explanation can be found in Thompson (1971, 1980) who proposed that this Rossby 
wave radiation mechanism is responsible for the emergence of strong eastward currents in the 
oceans and also conducted a laboratory experiment demonstrating the process (McEwan et al., 
1980). 
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Figure 1.17: Zonal jet emergence under localized statistically stationary forcing. The fluid 
has initially no mean flow. The red curve shows the linear jet amplitude growth that results 
in the mean from momentum convergence into the excitation region according to the classical 
wave-mean flow interaction theory, in which the modification of the eddy structure by the 
emerging mean flow is neglected. The blue curve shows the ensemble mean jet amplitude 
predicted to emerge by the statistical dynamical theory described in this thesis. The theory 
predicts exponential jet growth induced by the active eddy-mean flow interaction (dashed), 
followed by equilibration of the instability. The stochastic forcing consists of an ensemble of 
temporally delta-correlated waves with zonal wavenumber k x — 8 with Gaussian structure in y, 
i.e., e y ' , with d = 0.8. Other parameters are /5 = 1.4 and r = 0.1. 

The above theory assumes that the dominant and most relevant mechanism 
for jet emergence is the momentum convergence in the stirring region by the 
propagating Rossby waves. It is assumed that as the jet emerges its influence on 
turbulence (which is not all waves) is negligible and as a result it can, at first, 
be neglected (of second order if the jet is infinitesimal). However, we will show 
in this thesis and demonstrate immediately with an example, that the influence 
of the emerging jet on turbulence (which is neglected in the above theory) is 
the important and dominant mechanism for the emergence and maintenance of 
jets. Actually, it is dominant even for jets of infinitesimal amplitude. This active 
feedback of the mean flow on the turbulence results in a new type of instability 
that leads to exponential growth of the amplitude of the jet with the result that 
the amplitudes diverge exponentially from the linear growth predicted by classical 
wave-mean flow theory. It is important to note that this instability is an instability 
of the statistical dynamics of the turbulent flow. We will present in this thesis 
a second-order cumulant approximation to the full statistical dynamics of the 
turbulent flow that reveals this instability of the interaction between large-scale 
structure and turbulence. To demonstrate the implications of the statistical 
dynamical formulation of the wave-mean flow dynamics we plot in Fig. 1.17 the 
jet amplitude evolution predicted by the second-order closure theory discussed 
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in this thesis under the same forcing. The jet grows initially exponentially and 
then equilibrates, indicating that the second-order closure incorporates also the 
dynamics of equilibration. The instability manifests only in the ideal dynamics 
of the statistical state of the flow and is only partially reflected in individual 
simulations. For example, in Fig. 1.18a we show the development of the jet in a 
sample integration of the nonlinear equations of motion under a realization of the 
excitation that was used in Fig. 1.17 and in Fig. 1.18b a snapshot of the vorticity 
field and of the latitudinal structure of the jet that emerges. The nonlinear 
simulation confirms that the growth is faster than linear at first, but this sample 
integration can neither establish that there is an underlying instability nor make 
analytic predictions as what jet structure is expected to emerge at first or the 
parameter range that leads to jet emergence. Although this instability is revealed 
only when the dynamics of the statistical state of the turbulent flow are examined, 
the predictions of the statistical theory are reflected in sample simulations of 
the nonlinear dynamics. Moreover, this instability of interaction that leads to 
the emergence of jets does not even require that the forcing be localized. Jets 
may emerge even if the forcing is homogeneous, contrary to the predictions 
of classical wave-mean flow theory. In Fig. 1.19 we demonstrate in sample 
nonlinear simulations the emergence of robust jet structure both under spatially 
inhomogeneous forcing (Fig. 1.19a) and, most importantly, under homogeneous 
forcing (Fig. 1.19b). 

In this thesis we will use a non-equilibrium statistical theory to address for¬ 
mation and maintenance of jets and large-scale structures in turbulence. The 
proposed theory differs greatly from current theories that involve turbulent cas¬ 
cades and it has its basis in wave-mean flow interaction theories which consider 
that the most important interaction is the non-local in wavenumber space in¬ 
teraction between large-scale flows and the smaller scale eddies. Systematic 
investigation of the energy and enstrophy transfers among spectral components 
in numerical simulations has revealed that indeed the upgradient energy transfer 
from the small scales to the large-scale flow is mainly due to the highly non-local 
interactions in wavenumber space with a clear scale separation between them 
(Shepherd, 1987; Huang and Robinson, 1998). In this thesis we will demonstrate 
that not only local wavenumber interactions are not the main contributors to 
large-scale structure formation but moreover, they are not even required. 
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Figure 1.18: Reflection of the statistical dynamical results shown in Fig. 1.17 obtained by 
integrating the nonlinear barotropic equation (1.6) under a single realization of the forcing used 
in Fig. 1.17. (a): the time development of the jet amplitude, (b): a snapshot of the vorticity field. 
The realization reflects the predictions of the statistical theory that initially the jet amplitude 
grows exponentially. 



Figure 1.19: (a) Instantaneous snapshot of the vorticity field resulting from localized stirring 

in the region centered at y — 7r, indicated on the left of the panel, together with a snapshot of the 
structure of the zonal mean velocity (thick white line). This case differs from that of Fig. 1.18 in 
the spectrum of the excitation. Here the excitation is obtained by convolving a homogeneous 
excitation with isotropic spectrum centered at total wavenumber kf — 15 with a Gaussian 
localizing the excitation only in the y direction to the region around y — n. The dependence 
of the large-scale flows that emerge on the spectrum of the excitation will be addressed in this 
thesis. The jets in this case form from the momentum convergence into the excitation region 
resulting from the Rossby wave propagation and from the active feedback of the jet on the waves 
that leads to the intensification of the preexisting jet. (b) Snapshot of the zonal mean flow and 
of the vorticity field when the excitation is spatially homogeneous. While in this case classical 
wave-mean flow theory predicts that no mean flow should emerge, the statistical dynamics 
discussed in this thesis predict that a bifurcation occurs as the energy input, £/, increases or the 
dissipation coefficient, r, decreases. In the specific case, for parameter values £fk 2 /r 3 < 3.3 x 10 3 
the turbulent flow remains homogeneous with no jets, while for £fk 2 /r 3 > 3.3 x 10 3 the turbulent 
flow transitions to an inhomogeneous state with large-scale jets. The nonlinear integration shown 
here at £f k 2 /r 3 = 5 x 10 4 demonstrates that the predictions of the statistical dynamics are 
reflected in individual realizations of the flow. It will be demonstrated in this thesis that the 
predictions of the statistical dynamical theory are reflected also for parameters near the critical. 
Other parameters: /3/(kfr ) = 70 and the numerical integration was performed at resolution 
512 x 512 with a pseudospectral code. 28 
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Formulation of the S3T statistical state 
dynamics of turbulent flows on a /3-plane 


The formation and maintenance of zonal jets in planetary atmospheres is es¬ 
sentially governed by barotropic processes. The simplest setting in which we 
can study planetary barotropic processes is a planar flow on a rotating /3-plane 
which conserves the absolute vorticity of the flow in the absence of dissipation. 
Turbulence on a /3-plane does not self-sustain and a turbulent state must be 
externally forced in order to be maintained against dissipation. This forcing 
may model processes absent from the 2D barotropic dynamics, such as energy 
injected by baroclinic instabilities or turbulent convection. Because of the erratic 
and unpredictable nature of these vorticity sources in planetary turbulence, the 
forcing is modeled as a white-noise process in time given that the fluctuations of 
the forcing have a short autocorrelation time compared to the time scales of the 
barotropic dynamics. We also assume that the forcing is spatially homogeneous 
and if the turbulent flow becomes inhomogeneous this should be attributed to 
the dynamics. 

In the following sections we formulate the quasi-linear approximation of the 
nonlinear stochastically forced barotropic vorticity equations and derive the 
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equations of the S3T statistical dynamics of the turbulent flow on a barotropic 
/3-plane. S3T is an acronym for Stochastic Structural Stability Theory, which was 
initially abbreviated as SSST. The reason for this acronym will become apparent 
in this chapter. 

2.1 Formulation of the S3T dynamics on a /3-plane 

Consider a non-divergent, barotropic flow on a infinite /3-plane with planetary 
vorticity gradient, (3 — (0,/3). The velocity field being non-divergent can be 
expressed in terms of a streamfunction, as u = z x which implies (u^v) = 
(—dy^ip^dx^) (z is the unit vector normal to the /3-plane, see Fig. 1.10). The 
vorticity of the fluid is V x u = ( z with ( d x v — d y u — A i/j, and A = V • V = 
®xx + dy y the two-dimensional Laplacian. In the presence of stochastic forcing 
and dissipation, the potential vorticity, q = £ + /q + (3 • x, which here is simply 
the absolute vorticity, evolves as: 


— = d t ( + J(ip,C + P • x) = -r( + , 


(2.1) 


where D/D t = d t + u • V is the material derivative along the fluid flow. The 
advection term, (u • V)g, is alternatively expressed as where J(g,h) = 

(d x g)(d y h) — ( d y g)(d x h ) is the Jacobian of functions g and h. The flow is dissipated 
with linear damping at a rate r, which typically models Ekman drag in planetary 
atmospheres. Turbulence is maintained by the external stochastic forcing t ). 

We assume that is a homogeneous random stirring and we model this 
excitation as temporally delta-correlated Gaussian process with zero mean, i.e., 
(£(x, £)) = 0, and with spatial correlation prescribed by Q, 



( 2 . 2 ) 


The brackets denote the ensemble average over forcing realizations. (For details 
regarding the stochastic excitation refer to Appendix A.) We render (2.1) non- 
dimensional using as a time scale the dissipation time scale 1 jr and as a length 
scale the typical length scale of the stochastic excitation, Lf. With this non- 
dimensionalization (2.1) becomes an equation for the variables 



(2.3a) 
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and parameters 

' j *= 4 t ' e ‘ =;3 %' r * =1 - <2 - 3b) 

From here on we will work with the non-dimensional equation and drop the 
asterisks. Typical values for the non-dimensional parameters are /?* = 15, <s* = 
1200 for the Earth’s atmosphere, /?* = 3, £* = 1600 for the ocean and /?* = 450, 
£*=4 x 10 7 for Jupiter, based on the parameter values of the table 2.1. 

Table 2.1: Typical parameter values for geophysical flows. The typical forcing length scale is 
taken as the deformation radius in each geophysical setting. 



1 /kf 
[km] 

l/r 

[day(= 24h)] 

^rms 

[ms -1 ] 

/3 

[10 1:L m 1 s 1 ] 

fm- 2 V 3 ] 

/3* 

£* 

Earth’s atmosphere 

1000 

10 

15 

1.6 

2 x 10 -3 

15 

1300 

Earth’s ocean 

20 

100 

0.1 

1.6 

10 -9 

3 

1600 

Jovian atmosphere 

1000 

1500 

50 

0.35 

0.5 x 10 -5 

450 

4 x 10 7 


The first step in constructing the S3T dynamical system is to decompose the 
vorticity flow field into an averaged or mean field, Z — T[ C ], and deviations from 
the mean vorticity, (' = ( — Z, which is referred to as eddy vorticity. The averaging 
operator T determines the type of mean field we want to study. We employ 
two types of averaging operators: i) an average over the zonal x direction, i.e., 
T[(/)] = L~ x J^ x dx 0(x, £) and ii) a Reynolds average in which T[(j)\ produces a 
coarse-grained field which is obtained by averaging over an intermediate time scale 
or length scale which is larger than the time scale or length scale of the turbulent 
motions but smaller than the time scale or length scale of the coarse-grained field. 
In the first interpretation of the averaging operator the mean flows are zonal 
jets while in the second they may be either zonal jets or slowly moving traveling 
waves. 

With this decomposition, the barotropic vorticity equation (2.1) is equivalently 
rewritten as a system for the joint evolution of the mean and the eddy vorticity: 

dtZ + J0F,Z + /3-x) = -T[J(^,C')] -Z , (2.4a) 

dtC = -J (t//, Z + (3 • x) - J (* ,CO - C' + T[ J (<//, C)]-J (</>', CO + , 

'-V- ' V - ^' 

-4(U) C' /nl 

(2.4b) 

where T = T [ ] is the mean streamfunction and ?// is the eddy streamfunction. 

Equations (2.4) are referred to as the NL system. The stochastic excitation is 
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assumed to have zero mean, T[£] =0, and consequently the mean equations are 
unforced. The first term on the r.h.s. of (2.4b), w4(U)C', represents advection of 
eddy vorticity by the mean flow and is a bilinear functional of the eddies and 


the mean flow, while the second nonlinear term, f n \ = T [ J ('0 / , (') ] — J ('0 / , CO 


represents advection of the eddy field by itself. The operator *4(U), which governs 
the linear dynamics of the eddy field if the mean flow U = z x V$ is prescribed, 
can be written as: 



(2.5) 


We also make the ergodic assumption that the T average of a flow field 
(i.e. the zonal average or the Reynolds average over the intermediate time or 
length scale) is equal to the ensemble average over the forcing realizations, i.e., 
T[0(x, £)] = (</>(x,£)), where the brackets denote the ensemble average. The 
identification of the ensemble average with an averaging operation is crucial for 
the realization of the statistical quantities in a single planet and the validity of 
the ergodic assumption is established by experiment. 

In order to obtain a closed statistical description of the turbulent flow we 
restrict the nonlinearity in the NL equations by neglecting the eddy-eddy term, 
/ n i, in (2.4b) or parametrize it as stochastic excitation. We obtain in this way 
the quasi-linear (QL) approximation to the NL system (2.4): 


d t Z + J(T,Z + /3-x) = -r[J(^ / ,C / )] -Z , 
dtC = A(U) C' + ^C- 


(2.6a) 

(2.6b) 


A schematic comparing the nonlinear interactions operating in NL and QL system 
is shown in Fig. 2.1. In NL the term/ n i, which is neglected in QL, neither injects 
nor dissipates energy (see Appendix A) and therefore the QL system, in the 
absence of forcing and dissipation, has the same invariants as the NL system, 
namely it conserves both energy and enstrophy. 

The QL system has the attribute that its statistical dynamics close at second 
order. To obtain the statistical dynamics of the quasi-linear system (2.6) we 
use the ergodic assumption to identify Z — (() and the second cumulant of the 
vorticity between points x a and x^, 



(2.7) 
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Figure 2.1: Schematic of all possible wavenumber triad interactions in the NL and QL 
systems, (a) Two eddies with wavenumber vectors ki and k 2 combine to form a mean flow with 
wavenumber ki + k 2 . This interaction is in both NL and QL (it is term T[J ('ip',C) ] in (2.4a) 
and (2.6a)). (b) A wavenumber ki mean flow interacts with a k 2 mean flow to produce a mean 
flow with wavenumber ki + k 2 . This interaction is in both NL and QL (it is term J(T,Z) 
in (2.4a) and (2.6a)). (c) A wavenumber ki mean flow interacts with a wavenumber k 2 eddy to 
produce a ki + k 2 eddy. This interaction is in both NL and QL (it is term A(U)C / in (2.4b) 
and (2.6b)). (d) An wavenumber ki eddy interacts with a k 2 eddy to produce a ki + k 2 eddy. 
This interaction is included in NL (term f n \ in (2.4b)) but neglected in QL. 


with TlCfcait) C'( x fr,£) ]. Then, the average T[ CO] — can be 

expressed as a linear functional of C a b(t ) = C(x a ,Xb,£). To show that we use the 
incompressibility condition to rewrite J(?//, (') = V • [(z x V^') ('], and proceed 
as follows: 


r{v-[(zxv^)ci} 


v-r [(z x vV’Oc'] = v • ((z x O 


V-(-ZX (V 0 <Cb + V 6 ^Ca) 


zxfvAbv^ 1 ) 


x«=x 6 

C ab 


Xa=*b 


(2.8) 


The subscript a or b in functions denotes hereafter the value of the function at 
x a or x&, i.e. (J a = £'(x a ,£), the subscript a or b in operators denotes the action 
of the operators only on the variables x a or x^ respectively, and the subscript 
x a = X6 denotes that any expression depending on the two variables x a and 
x^ should be evaluated at x a = x^ = x. 1 The operator A -1 is the inverse of 
the Laplacian which has been rendered unique by incorporating the boundary 
conditions. Equation (2.8) shows that T[J('0 / , CO] is a linear functional of C. We 
denote the linear functional given in (2.8) by 1Z and set 


K{C) = -T[J( YX')] • (2-9) 

J For example: d y {d Xb [C'(x 0 , t)C'(x b , i)]} = d y [C'(x a , t) d Xb C(*b, *)] = 

V. J X a —Xf, L J x a —Xf, 

dy X [C'tM) 2 / 2 ]- 
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The first cumulant, Z, of the flow field therefore evolves according to: 

d t Z + J O, Z + (3 • x) = K{C) - Z . (2.10) 

To obtain the evolution equation of C we take the time derivative of (2.7) and 
obtain: 2 


d t C ab = [A(U) + A(U)] C ab + \fe (Ub + CUb 


( 2 . 11 ) 


where *A a (U) indicates that the coefficients of the operator *4(U) are evaluated 
at x a and that the differential operator act only on the variable x a of C(x a ,Xb,£). 
(Similarly for 7U(U).) 

It can be shown (see Appendix A. 2) that for temporally delta-correlated 
stochastic forcing term (^aCb + CaCfr) is independent of the state of the system 
and exactly equal to — x&) = y/eQab 3 The rate of energy injection 

is thus independent of the state of the system and is prescribed by the spatial 
forcing covariance Q and the amplitude factor e. The same is true for the NL 
system. In both systems the energy injection rate is e 
where Q( k) is the Fourier transform of Q, 


(2tt) -2 / d 2 k Q(k)/ (2/c 2 )J. 


Q(k) = j 


d 2 (x a -x ft ) Q(x a -x b )e 


-ik-(x a -x b ) 


( 2 . 12 ) 


with k = (k x ,ky). Because Q is a homogeneous covariance its Fourier transform is 
real and non-negative, i.e., Q( k) >0, for all wavenumbers k. The quantity Q/k 2 , 
where k = |k|, determines the energy spectrum of the forcing. We normalize Q so 
that 

f d 2 k Q(k) = 

J (2t r) 2 2 k 2 

and the energy injection rate per unit area is e. For details see Appendix A. 

The joint evolution of the first two cumulants of the flow field, Z and (7, define 
the S3T statistical state dynamics of the turbulent flow which is governed by the 

2 In writing (2.11) we adopt the Stratonovich interpretation for stochastic differential equations. 
However, because the stochastic forcing in our case is additive, both Stratonovich and Ito 
interpretations lead to the exact same results (cf. Appendix A). 

3 The dependence of the spatial covariance of the forcing on the difference coordinate x a — 
indicates that the stochastic forcing is spatially homogeneous (see Appendix A). 
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autonomous system of deterministic equations: 


(2.13a) 

(2.13b) 


d t Z + J (*, Z + (3 ■ x) = K{C) - Z , 


dtCab 


A a (U) + A b (U)]c ab + £Q ab . 


The S3T system (2.13) corresponds to a second-order closure of the full sta¬ 
tistical dynamics of the turbulent flow. This closure became possible because of 
the adoption of the quasi-linear approximation. If the quasi-linear approximation 
were not made, then the evolution of the second cumulant, C, would also involve 
terms of the form (/ n i(x a , t)(' / (x 5 , £)), which are related to the third cumulant and 
as a result the equations for the first two cumulants would not close. Neglecting 
or parametrizing the eddy-eddy terms in (2.4b) by a state independent Gaussian 
stochastic process leads to a closed set of equations for the evolution of the first 
two cumulants of a Gaussian approximation of the statistical state dynamics of 
the turbulent flow. This approximation is also referred to as “CE2”. Note that 
higher order truncations of the cumulant equations is problematic. Marcinkiewicz 
(1939) has shown that truncations of the cumulant equations at order no > 2, 
obtained by setting all n-th order cumulants for n > no equal to zero, produce 
non positive probability density functions (pdf). Therefore the only physically 
realizable cumulant truncation is at second order. 

2.2 Formulation of the S3T dynamics of zonal mean states 

The most common mean flows that appear in planetary turbulence are zonal 
jets. In order to address the statistical dynamics of zonal jets in turbulence we 
may choose the averaging operator T to be directly the average over the zonal 
direction, x, i.e., 

T[(f>} = j-[ (j>(x',y,t)dx' = $(y,t) . (2.14) 

L x J o 

The zonal average of a flow field is also denoted with an overbar, for example 
0(x, t) = L x ~ x Jq x (j){x f ,y,t) dx f . This zonal S3T closure simplifies significantly 
the QL and S3T systems and it is the easiest to interpret because the separation 
between mean and eddy is unequivocal. 

With the zonal average the mean flow vorticity Z is related to the zonal flow 
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mean flow, £/, through Z — —d y U, and because non-divergence implies d y V = 0, 
without any loss of generality, we can assume that V = 0. Since /3 — (0, /3) the 
advection of the mean potential vorticity flow, Z + /3 • x, by the mean flow field, U, 
vanishes, i.e., J (T, Z + {3 • x) = 0, and the zonal average vorticity flux divergence 


simplifies to: 

J (V, CO = v • (u' CO = d y (yc) • (2.15) 

With these simplifications the NL system takes the form: 

d t u = VC -U , (2.16a) 

dtC = A Z (U) C' + d y (VC 7 ) - V • (u'CO , (2.16b) 

fnl,z 

while the QL system becomes 

d t U = ¥C~U , (2.17a) 

dtC = A Z {U)C + V~£^ (2.17b) 

where in both (2.16) and (2.17) operator A z is 

MU) = -Ud x -(p- d 2 yy U ) dxA- 1 - 1 . (2.18) 


(Roman subscript z denotes that the zonal mean-eddy decomposition was used.) 

The three types of nonlinear triad interactions that can occur between the 
mean quantities and the eddies are shown in Fig. 2.2. In the QL approximation we 
neglect / n i ?z or parameterize it as stochastic noise. It should be noted that Bouchet 
et al. (2013) have established that in this zonal mean-eddy decomposition the 
QL approximation becomes exact in the limit of £* = £/(r 3 Lj) —oo (cf. (2.3b)). 

Because A Z (U) is invariant under the translation x —>> x + a for any constant <a, 
the eddy vorticity equation is homogeneous in x and therefore the eddy vorticity 
covariance will always be homogeneous in x and consequently of the form: 


C(x n , x 6 , t) = C(x a - x b , y a , y b , t ) . 


(2.19) 


The zonal homogeneity of C allows us to simplify the flux divergence to: 


7 Z Z (C) = -dy 


(^a d Xa +A b d Xb ^j C ab 


( 2 . 20 ) 


X a =X(, 
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Figure 2.2: Schematic of the triad interactions among the zonal or x -wavenumber fields in 
the NL system (2.4). The wavenumbers in this figure refer only to the zonal or x-wavenumber 
components of the full wavenumber vector k. (a) Two eddies with wavenumbers k x and — k x 
interact to form a zonal mean flow, U (with k x = 0). This interaction is in both NL and QL 
(it is term v'(' in (2.16a) and (2.17a)). (b) A k x wavenumber eddy interacts with mean flow U 
(k x = 0) to produce another eddy with zonal wavenumber k x . This interaction is in both NL 
and QL (it is term A Z (U) (' in (2.16b) and (2.17b)). (c) A k x \ wavenumber interacts with a k x 2 
wavenumber eddy to produce a k x \ + k x 2 wavenumber eddy. This interaction is included in NL 
(term / n i, z in (2.16b)) but neglected in QL. 


and the S3T system takes the form 

1 


d t U = 


' A X v Xh )^ab 

Xa=*b 

dfCab = A z ,a(U) + A z , b (U) ^ab + £ Qab • 


-U, 


(2.21a) 

(2.21b) 


This system will be denoted as S3Tz (for S3T-zonal). 

Solutions of (2.21) are also solutions of the generalized S3T system (2.13), i.e., 
a solution (u(y, £), C(x a — y a , t)^ that satisfies (2.21) also satisfies (2.13) 
as well; the converse however is not true. S3Tz system (2.21) has tremendous 
advantage in numerical simulations over the generalized S3T system (2.13) because 
its state variables have significantly fewer degrees of freedom. The method of 
numerical integration of the stochastic NL and QL and of the deterministic S3T 
equations is discussed in Appendix C. 


2.3 S3T STATISTICAL EQUILIBRIA AND THEIR STABILITY 

S3T systems (2.13) and (2.21) are autonomous and may admit equilibrium (fixed 
point) solutions ^Z e (x), C e (x a , . These equilibria are statistical equilibria of 
the turbulent flow and consist of the mean flow vorticity Z e {x ) and an eddy field 
with covariance C e (x a ,xj ) ). 

Remarkably, both S3T systems (2.13) and (2.21) admit the stationary homoge- 
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neous equilibrium 


( 2 . 22 ) 


Z e = 0, C e = |Q, 

for all values of e and /3 — (/3 x ,/3 y ) with components /3 X and /? y , under the 
condition that the forcing covariance is homogeneous. This statistical equilibrium 
has no mean flow and a homogeneous eddy field. To confirm this note that 
= A{ U e = 0) = -1 + z • (/3 x V)A _1 , and 


(A e a + At)C e = { - 2 + z • [/3 x (V.A- 1 + V^A," 1 )] } e -Q 

\f3x (VaA- 1 + V b A^)] I d ' k 


= ~eQ + 2% ' 


Q(k)e ik ' (x “- X6) 


= -eQ - 
= -eQ , 


/ 


d 2 k 
(27r) 2 


(3 


ik -ik 

+ 


— k 2 —k 2 


(27r) 2 
Q(k)e ik '( x “- x ^ 


(2.23) 


showing that C e of (2.22) satisfies (2.13b). Further, from (2.9), 

d 2 k 


7Z(C e ) = --V■ 


Z X 


Z X 


(v.A- 1 + v&a 6 

d 2 k ( ik -ik 


'A / (2rf 


/ 


+ 


(27r) 2 \—k 2 —k 2 


(27T) 2 
Q(k)e ik -( Xa " x ^ 


X a =X 6 

= 0, 


J x a =x b 


(2.24) 


which in turn confirms that (2.13a) is also satisfied. While the homogeneous 
state (2.22) is always an equilibrium of the S3T system it may only be an 
approximate equilibrium of the full hierarchy of cumulant equations. However, 
we show in Appendix G that for the case of isotropic delta function ring forcing, 
i.e., for Q( k) = Ank 2 8{k — kf), this homogeneous statistical equilibrium is also 
an equilibrium of the full hierarchy of cumulants. 

The stability of any S3T equilibrium solution (Z e ,G e ) is addressed by con¬ 
sidering small perturbations (SZ, SC) about this equilibrium and performing an 
eigenanalysis of the linearized S3T equations about this equilibrium: 

d t 5Z = A e SZ + K(6C) , (2.25a) 

d t SC ab = (A e a + At) 8C ab + (« 5Aa + SAi,) C e ab , (2.25b) 

where A e = A( U e ) and 5A = A{ U e + <fU) - A e . 
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When the equilibrium is unstable the statistics of the flow bifurcate to a new 
state. So the stability of an S3T equilibrium implies the structural stability of 
the turbulent flow, while the marginally stable S3T equilibria identify the critical 
states at which the turbulent flow becomes structurally unstable and transitions 
to a new statistical state. 

S3T stability involves the stability of the statistics of the turbulent state and 
is fundamentally different from the hydrodynamic stability of a mean state. It 
can be shown that if the S3T equations admit the equilibrium ( Z e ,C e ) then by 
necessity the associated mean state is hydrodynamically stable (cf. Appendix B). 
However, the hydrodynamic stability of a mean state does not imply the S3T 
stability. Most notable example is the homogeneous equilibrium with no mean 
flow (2.22). The state of zero mean flow is clearly hydrodynamically stable but it 
will be shown that at a critical parameter e the homogeneous equilibrium becomes 
S3T unstable and the turbulent flow reorganizes to an inhomogeneous state. 

2.4 Bibliographical Note 

The S3T theory was introduced by Farrell and Ioannou (2003). The continuous 
formulation of the theory was developed by Srinivasan and Young (2012). The 
cumulant interpretation was discussed by Marston et al. (2008) who refer to it as 
CE2 (see also Marston (2012)). The cumulant representation of the statistical 
dynamics of the flow were developed by Hopf (1952). The statistical stability of 
the homogeneous state in S3T (or CE2) and the subsequent formation of zonal 
jets is investigated in barotropic flows by Farrell and Ioannou (2007); Bakas 
and Ioannou (2011); Srinivasan and Young (2012); Parker and Krommes (2014a). 
Earlier, Carnevale and Martin (1982) using field theoretic techniques arrived at the 
same equations for the statistical description of fluctuations about a homogeneous 
state but the relevance for the emergence of zonal jets was not discussed. The 
statistical stability of inhomogeneous states in S3T is investigated by Farrell and 
Ioannou (2003); Parker and Krommes (2014a). Statistical state dynamics with 
higher order cumulant truncations are discussed by Marston (2012); Marston 
et al. (2014). The generalized coarse-grained mean flow interpretation of S3T that 
allows non-zonal solutions was introduced by Bernstein and Farrell (2010) in an 
investigation of the phenomenon of blocking in a two-layer baroclinic atmosphere 
and was studied recently for barotropic flows by Bakas and Ioannou (2013a, 2014). 
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3 

Emergence of coherent structures out of 
homogeneous turbulence through S3T 

instability 


3.1 S3T INSTABILITY OF HOMOGENEOUS TURBULENT EQUILIBRIUM 

We have seen in the previous chapter that for spatially homogeneous forcing there 
is always a homogeneous equilibrium of the S3T system (2.13). This equilibrium 
is given by 

Z e = 0 , C e (x a - x 6 ) = |<2(x a - x fc ) . (3.1) 

We want to determine the statistical stability of this equilibrium as a function 
of the parameters available in the problem. These parameters are the non- 
dimensional e and /? defined in (2.3b) and also the spectrum of Q. We examine 
cases in which the spectrum of Q is isotropic and cases in which it is anisotropic. 
The stability of this equilibrium is determined by the linearized S3T perturbation 
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equations (2.25) about the homogeneous equilibrium (3.1), which take the form: 


d t SZ = A e 5Z + K(i5C) , (3.2a) 

dt 8C a b = (Al + At) 8C ab + (SAa + 8A b ) C e ab , (3.2b) 


where (SZ, SC) are the perturbation mean flow and perturbation covariance, 
A e = z • ({3 x V) A" 1 - 1 and SA = -6U • V + f(A<SU) • vl A 


-l 


The purpose of this chapter is to examine the stability of the homogeneous 
equilibrium. We derive an analytic expression for the eigenvalues of (3.2) and 
show that there is always a critical energy input rate e = £ c that renders (3.1) 
unstable. When the equilibrium is unstable a mean flow in the form of the most 
unstable mean flow eigenfunction grows, initially at the rate predicted by the 
eigenvalue, and the turbulent flow will eventually reorganize to an inhomogeneous 
state. We study the dependance of £ c on non-dimensional /? for isotropic and 
anisotropic forcing spectra and also determine which type of mean flow (zonal 
jets or non-zonal flows) is the most unstable. 


3.2 Eigenanalysis of the homogeneous equilibrium 

We proceed now with the stability analysis of (3.1). Consider eigenfunctions of 
the form (SZ,SC) e st . The eigenvalue s and 5Z(x) and 5C(x a ,X6) satisfy the 


eigenvalue problem 

s5Z = A e SZ + n(5C ) , (3.3a) 

s 5C ab = (At + At) 5C ab + (8Aa + 8A b ) C e ab . (3.3b) 

The eigenfunctions can be assumed in the form 

5Z n (x) = e in ' x , (3.4a) 

5C n (x a , x 6 ) = (7„ h) (x a - x 6 ) e in '( x <x+x (,)/2 } (3.4b) 


with n = ( n x ,n y ) the wavevector of the eigenfunction (see Appendix E). The 
mean flow component of the eigenfunction (3.4a) is a zonal jet when n x — 0 and 
a non-zonal flow, a plane wave, when n x A 0. 

Note that the mean flow eigenfunction SZ n is also an eigenfunction of A e with 
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eigenvalue — (ia; n + 1 ) 


(3.5) 


A e 8Z = -(kj n + l)6Z , 
where oo n is the Rossby frequency 




z ■ (18 x n) 


n* 


(3.6) 


with n — |n| and as a result (3.3a) can be written as (cr + 1)8 Z — 7 Z(8C) with 
a = s + icj n . Because the Reynolds stress associated with the perturbation 
covariance, 7 Z(8C), is proportional to e (as C e is proportional to e) it can be 
written as 

7 Z(8C) = ef(a)8Z , (3.7) 

where / is the Reynolds stress feedback or eddy feedback, and cr satisfies the 
dispersion relation 


cr + 1 =ef(a) 


(3.8) 


We remind the reader that the 1 in the l.h.s. of (3.8) is the rate of dissipation and 
therefore the homogeneous state is unstable when Re(cr) > 0 or £:Re /(cr) > 1, 
i.e., the mean flow acceleration by the Reynolds stress feedback exceeds the decay 
due to dissipation. The term f(a) measures the feedback on the mean flow 8Z 
by the eddy perturbation field after being distorted by the mean flow 8Z. When 
Re [/(cr)] > 0 the feedback on the mean flow by the eddy perturbation field has 
the tendency to reinforce the existing mean flow and the vorticity fluxes due to the 
eddies are upgradient. It is necessary for instability to have upgradient vorticity 
fluxes but it is not sufficient, because they have to overcome the dissipation. In 
Appendix E we show that the function f(cr) is (cf. Appendix E, eq. (E.10)): 


IW) = / 


d 2 k 

(2^F 


|k x n | 2 (fc 2 — k 2 )(k 2 — n 2 ) 


0(k) 


k^k^n 2 


(cr + 2) + i (cJk+n — ~ ^k) 


(3.9) 


with k s = k + n and k s = |k s | and the dispersion relation for the stability of the 
homogeneous equilibrium is 


cr + 1 


= £ J 


d 2 k 

( 2 vr ) 2 


|k x n | 2 (kg — k 2 )(k 2 — n 2 ) 


Q(k) 


k^k^n 2 


(<7 + 2) + i (Wk+ n — iO n — CJk) 


(3.10) 
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Figure 3.1: Top panels: the forcing covariance spectrum, Q(k) = An 5(k — 1) [1 + fj, cos(2j)], 
for (a) fi = 1, (b) /i — 0 and (c) /i = — 1 (the support of the delta function is represented as a thin 
ring). Bottom panels: contours of the vorticity field induced by a realization of the stochastic 
forcing for (d) fi = 1 , (e) [i = 0 and (f) n — — 1 . 

We investigate the stability of the homogeneous equilibrium under stochastic 
forcing with spectrum 

Q(k) = 47r ) 5(k-l) , (3.11) 

with 7 = arctan (k y /k x ) and 

Q( 7 ) = 1 + /i cos(27) . (3.12) 

This forcing excites an eddy field at total wavenumber kf — 1 (in dimensional 
units kf = 1/Lf) and the parameter /i measures the anisotropy of the forcing. 
Parameter (i takes values \fi\ < 1 so that ^( 7 ) > 0 for all 7 . For /x = 0 the forcing 
is isotropic (see Fig. 3.1b,e). For /x > 0 the stochastic forcing is anisotropic (see 
Fig. 3.1a) favoring structures aligned with the meridional axis (i.e. with k y — 0), 
as shown in Fig. 3.Id, while /1 < 0 (see Fig. 3.1c) favors structures aligned with 
the zonal axis (i.e. with k x = 0), as shown in Fig. 3.If. In Jupiter because the 
excitation models vorticity input by turbulent convection we expect excitation to 
be of the fi — 0 type, while in the Earth because the excitation models injection 
of vorticity due to baroclinic processes we expect excitation closer to /i = 1 . 
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We determine the critical energy input rate £ C5Z that renders the homogeneous 
equilibrium unstable to zonal jet perturbations and the critical energy input rate 
£c,nz that renders the homogeneous equilibrium unstable to non-zonal pertur¬ 
bations. £ C;Z is the minimum e for which Re(cr) = 0 for an eigenfunction with 
wavevector n = (0, n y ) and £ c , nz is the minimum e for which Re(cr) = 0 for an eigen¬ 
function with wavevector n = (n x , n y ) and n x ^ 0. When e > min (e CjZ , £ C}nz ) = £c 
the homogeneous equilibrium is unstable and the structure that first emerges is 
zonal or non-zonal according to whether the minimum e is £ C;Z or £ c ,nz- 

The critical energy input rates, £ CjZ and £ c ^ nz as a function of (3 for isotropic 
forcing (/i = 0) is shown in Fig. 3.2a. For /? < 3.5 the structures that become 
first unstable are zonal jets [e c ^ z < £ c ,nz) and for supercritical energy input rates 
always zonal jets are more unstable than non-zonal perturbations. For /? > 3.5 
non-zonal structures become first unstable and for a range of energy input rates 
£c,nz < s < £ c ,z only them are unstable. For s > e c ^ z zonal jets become unstable 
but with less growth rates compared to non-zonal structures. For /? < 3.5 zonal jet 
eigenfunctions grow the most whereas for /? > 3.5 non-zonal structures grow the 
most. In the light shaded region only non-zonal coherent structures are unstable, 
while in the dark shaded region both zonal jets and non-zonal coherent structures 
are unstable. Growth rates, oy, as a function of the eigenfunction wavevector 
n — riy) for 4 different choices of [3 and s are shown in Figs. 3.2b-e. 

The £ C5 z and £ c?nz for both isotropic as well as anisotropic forcing are shown in 
Fig. 3.3. This figure shows that the homogeneous equilibrium becomes unstable 
for all values of (3 > 0. The homogeneous equilibrium becomes also unstable 
even for (3 = 0, unless the excitation is exactly isotropic (cf. Appendix E.l). 
This is an important result because it shows that the dynamics that lead to the 
initial emergence of large-scale structure does not require the presence of /?. We 
show in the next sections that for isotropic forcing both £ C5Z and £ CjIlz increase 
as /3~ 2 as /? 0, but for anisotropic forcing £ c = 32/|/i| + 0(f3 2 ) for small (3. 

The homogeneous equilibrium is rendered unstable with the least e in the range 
1 < /3 < 10. For /? > 4 the equilibrium becomes first unstable to non-zonal 
perturbations. As (3 increases the homogeneous equilibrium becomes more stable 
and larger e is required to destabilize it. It is shown that zonal jet emergence 
requires £ CjZ ~ f3 2 as (3 oo, which means that the effective feedback on the 
mean flow falls as Re [/(a)] ~ (3~ 2 as /3 oo, but for the emergence of non-zonal 
structure £ c ^ nz ~ /3 1 / 2 as f3 oo (Re[/(cr)] ~ /3 -1 ^ 2 ) because of the occurrence 
of fortuitous resonances that are explained in the next sections. The asymptotic 
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behavior of £ C5Z and £ CjIlz for large /? is independent of the forcing spectrum. 


3.3 Eddy-mean flow dynamics underlying the S3T instability of 

HOMOGENEOUS TURBULENT EQUILIBRIUM 

In order to analyze the dynamics underlying the S3T instability we study the 
behavior of Re(/) at the critical e at which the eigenfunction with wavevector 
n becomes neutral and set a r = 0 and <j{ = 0 in (3.9) and (3.10). 1 We denote 
the eddy feedback on the mean flow perturbation with wavenumber n in this 
approximation as f r = Re(/(0)). The eddy feedback for this delta function 
forcing (3.11) can be written as 

7r 

/r = yV(0, n)d0, (3.13) 

0 

where n) is the contribution to f r from the individual forcing components of 
Q corresponding to wavenumbers k and — k. For the narrow ring forcing (3.11) 
all forcing components have k — 1 and are only characterized by angle 0 , that 
is subtended measured from the lines of constant phase of the eigenfunction n 
(see Fig. 3.4). We also write n = (rising, ncostp) so that zonal jet eigenfunctions 
correspond to ip = 0°, while non-zonal eigenfunctions to <p ^ 0°. The angle 
7 = arctan(fc y /fc :r ) is given as 7 = 6 — p>. The relation between angles 0, cp and 7 
is shown in Fig. 3.4. We can isolate the dependence of this eddy feedback on /? 
by writing as J r (0, n) = F(Q, n) + F(180° + 0, n) with 

F (0,n) = p 2 + /3 2 °p 2 > ( 3 - 14 ) 

where, as shown in Appendix E.2, functions A/*, Vq and V 2 do not depend on (3. 
T measures the feedback on the mean flow from two monochromatic excitations 
with wavenumbers k and — k (see Fig. 3.4). We wish to determine the 6 that 
produce positive feedback to eigenfunction n and contribute to the instability of 

n. 

In the following sections we determine the contribution of the various waves to 
the eddy feedback and identify the angles 6 that produces the most significant 
contribution to this feedback. We also calculate the eddy feedback f r as a function 

x That cfi — 0 or equivalently Si = — for all wavevectors n at the stability boundary is an 
approximation but it can be shown that is a valid approximation. 
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Figure 3.2: (a) The critical energy input rates £ CjZ (solid) and s c ,nz (dashed) that render 

the homogeneous equilibrium unstable to zonal jet perturbations (n = (0, n y )) and non-zonal 
perturbations (n = (n x ,n y )) respectively, as a function of p for isotropic forcing covariance 
spectra (fi = 0, Fig. 3.1b). Also shown are the slopes P~ 2 , f3 2 and p 1 ^ 2 (dash-dot). Typical values 
of the Earth’s atmosphere and ocean and Jupiter’s atmosphere (found in table 2.1) are marked 
with stars. For p < 3.5 zonal jet eigenfunctions grow the most whereas for P > 3.5 non-zonal 
structures grow the most. In the light shaded region only non-zonal coherent structures are 
unstable, while in the dark shaded region both zonal jets and non-zonal coherent structures are 
unstable, (b)-(e) S3T growth rates, cr r , as a function of the eigenfunction wavevector n = (n x ,n y ) 
for the four cases marked in (a). The thick line corresponds to the a r = 0 contour. For (b), (d), 
(e) the contour interval is 0.15 while in (c) the contour interval is 0.5. The dashed line marks 
n — 1 and corresponds to the forcing scale. 
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Figure 3.3: The critical energy input rates s c ,z (solid) and e c ,nz (dashed) that render the 
homogeneous equilibrium unstable to zonal jet perturbations (n = ( 0 , 71 ^)) and non-zonal 
perturbations (n = (n x ,n y )) respectively as a function of p. Shown are the £ c , z and £ Cjnz for 
the three forcing covariance spectra seen in Fig. 3.1. Also shown are the slopes /3 ~ 2 , /3 2 and 
/ 3 1//2 (dotted line). For p > 4 the equilibrium hrst becomes unstable to non-zonal perturbations 
regardless of fi. For /jl m — 1 zonal jet perturbations are unstable only for p > 1.8. 



Figure 3.4: A non-zonal plane wave perturbation with wavevector n at an angle ip to the 
northward direction (the direction of (3) becomes a zonal perturbation when the coordinate frame 
is rotated clockwise by ip. Under this rotation the components of the wavevector k = (cosy, sin 7 ) 
are transformed to k = (cos 0 , sin 6), with 0 = 7 + p. F(n, 9) in (3.13) is the mean momentum 
flux convergence from plane wave perturbations that arise from excitations with wavevectors k 
and —k. 
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of the total mean flow wavenumber n for 0 < p < 90°. We limit our discussion to 
the emergence of mean flows with n < 1, i.e., with scale larger than the scale of 
the forcing. (Remember that all wavenumbers are non-dimensionalized with the 
forcing wavenumber kf .) In section 3.4 the analysis is mostly focused to isotropic 
forcing (Q = 1) while the effect of anisotropy is discussed in section 3.5. 


3.4 Eddy-mean flow dynamics leading to formation of zonal and 

NON-ZONAL STRUCTURES FOR ISOTROPIC FORCING 

3.4.1 Induced vorticity fluxes when /3 < 1 
We expand the integrand T of (3.13) in powers of (3: 


F = + + 0(/3 4 ) , 


(3.15) 


p=0 


. The leading order term, T\ o, is the contribution of each 


with = \ 

wave with wavevector k = (cos #, sin 6) to the eddy feedback in the absence 
of f3 and is shown in Fig. 3.5a. For (3 — 0, the dynamics are rotationally 
symmetric and for isotropic forcing f r is independent of p. Therefore all zonal 
and non-zonal eigenfunctions with the same total wavenumber, n, grow at the 
same rate. Upgradient fluxes (J^ > 0) to a mean flow with wavenumber n are 
induced by waves with phase lines inclined at angles satisfying 4 sin 2 6 < 1 + n 2 
(cf. Appendix F). This implies that all waves with \0\ < 30° necessarily produce 
upgradient vorticity fluxes to any mean flow with wavenumber n < 1 , while 
waves with 30° < \6\ < 45° produce upgradient fluxes for any mean flow with 
large enough wavenumber (cf. Fig. 3.5a). The eddy-mean flow dynamics was 
investigated in the limit of n <C 1 by Bakas and Ioannou (2013b). It was shown 
that the vorticity fluxes can be calculated from time averaging the fluxes over the 
life cycle of an ensemble of localized stochastically forced wavepackets initially 
located at different latitudes. For n <C 1, the wavepackets evolve in the region 
of their excitation under the influence of the infinitesimal local shear of 5U and 
are rapidly dissipated before they shear over. As a result, their effect on the 
mean flow is dictated by the instantaneous (with respect to the shear time scale) 
change in their momentum fluxes. Any pair of wavepackets having a central 
wavevector with phase lines forming angles \0\ < 30° with the y axis surrender 
instantaneously momentum to the mean flow and reinforce it, whereas pairs with 
\6\ > 30° gain instantaneously momentum from the mean flow and oppose jet 
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formation. Therefore, anisotropic forcing that injects significant power into Fourier 
components with \0\ < 30° (such as the forcing from baroclinic instability that 
primarily excites Fourier components with 0 = 0°) produces robustly upgradient 
fluxes that asymptotically behave anti-diffusively. That is, for a sinusoidal mean 
flow perturbation 5U = sin (ny) we have Jq Tq d 9 = kvP with k positive and 
proportional to the anisotropy factor \i (cf. Appendix F). 

For isotropic forcing the net vorticity flux produced by shearing of the pertur¬ 
bations vanishes, i.e., Jq J-'o d9 = 0, given that the upgradient fluxes produced 
by waves with \6\ < 30° exactly balance the downgradient fluxes produced by 
the waves with \9\ > 30°. However, a net vorticity flux feedback is produced and 
asymptotically behaves as a negative fourth-order hyperdiffusion with coefficient 
0(/3 2 ) for /? <C 1 (cf. (3.16) and Bakas and Ioannou (2013b)). In Appendix F it 
is shown that the feedback factor f r for isotropic forcing in the limit /3 <C 1 with 
f3/n <C 1 is: 


fr 


>n 


/? — 2 + cos(2 ip) + 0(P ) , 


(3.16) 


which is accurate even up to n ~ 1, as shown in Fig. 3.6. In order to understand 
the contribution of /? to the vorticity flux feedback, we plot J^/n 4 for a zonal 
(Fig. 3.5b) and a non-zonal perturbation (Fig. 3.5c) as a function of the mean flow 
wavenumber n and wave angle 9. We choose to scale T 2 by n 4 because in (3.16) 
f r increases as n 4 . Consider first the case of a zonal jet. It can be seen that at 
every point, T 2 has the opposite sign to implying that /? tempers both the 
upgradient (for roughly \9\ < 30°) and the downgradient (for \9\ > 30°) fluxes of 
Tq. However, in the sector \9\ > 30° the values of T 2 are much larger than in the 
sector \9\ < 30° and the net fluxes integrated over all angles are upgradient, as 
in (3.16) for the isotropic case. 

The asymptotic analysis of Bakas and Ioannou (2013b), which is formally 
valid for n <C 1, offers understanding of the dynamics that lead to the inequality 
^ 2^0 < 0 and to the positive net contribution of i.e., to /q T 2 d9 > 0. Any 
pair of wavepackets with wavevectors at angles \0\ > 30° instantaneously gain 
momentum from the mean flow as described above (i.e. Jjj < 0 for \9\ > 30°), but 
their group velocity is also increased (decreased) while propagating northward 
(southward). This occurs due to the fact that shearing changes their meridional 
wavenumber and consequently their group velocity. The instantaneous change 
in the momentum fluxes resulting from this speed up (slowing down) of the 
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Figure 3.5: (a) Contours of JT(0, ra) in a (0, n) polar plot (n radial and 0 azimuthal). This 

figure shows the magnitude and sign of the vorticity flux induced by waves with phase lines 
oriented at an angle 0 to the y axis in the presence of an infinitesimal mean flow perturbation 
of total wavenumber n when /? = 0. The contour interval is 3 x 10 -3 and note that To (0, n) is 
independent of <p. (b) Contours of the normalized Ti (0,n)/n 4 show the 0(/3 2 ) correction to 
JT(0, n) for the case of zonal jet perturbations (<p = 0°). The contour interval is 0.02. (c) Same 
as (b) but for non-zonal perturbations with <p — 15°. The contour interval is 0.04. In all panels 
the forcing is isotropic (fi = 0), solid (dashed) lines indicate contours with positive (negative) 
values, the thick line is the zero contour, the radial grid interval is An = 0.25 and the 30° wedge 
is marked (dashed-dot). In panels (a) and (b) the zero contour is the curve 4sin 2 0 = 1 + n 2 
(see Appendix F). 
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Figure 3.6: Eddy feedback factor f r as a function of n for /3 = 0.1 and isotropic forcing. Note 
that the fluxes are upgradient (i.e. f r > 0) for all mean flow wavenumbers n. Shown is f r for 
ip = 0° and p = 60° (solid lines), as well as the asymptotic expression (3.16) (dash-dot) derived 
for the feedback factor in the limit <C 1 and /3/n <C 1. 

wavepackets is positive in the region of excitation leading to upgradient fluxes 
(T 2 > 0). The opposite happens for pairs with \0\ < 30° (cf. Fig. 3 of Bakas 
and Ioannou (2013b)), however the downgradient fluxes produced are smaller 
than the upgradient fluxes, leading to a net positive contribution when integrated 
over all angles. Figure 3.5b, shows that this result is valid for larger mean flow 
wavenumbers as well. 

Consider now the case of a non-zonal perturbation (Fig. 3.5c). We observe that 
the angles for which the waves have significant positive or negative contributions 
to the vorticity flux feedback are roughly the same as in the case of zonal jets. 
In addition, the vorticity flux feedback factor decreases with the angle p of 
the non-zonal perturbations (cf. (3.16)). As a result, zonal jet perturbations 
always produce larger vorticity fluxes compared to non-zonal perturbations and 
are therefore the most unstable in the limit f3 1. Additionally, these results 
show that for /? <C 1, the mechanism for structural instability of the non-zonal 
structures is the same as the mechanism for zonal jet formation, which is shearing 
of the eddies by the infinitesimal mean flow. 

3.4.2 Induced vorticity fluxes when (3 » 1 

When /? >> 1 by inspecting (3.14) we expect that f r should fall as /3~ 2 . This indeed 
is the case, as we will demonstrate, for zonal jet perturbations. However, non-zonal 
perturbations may render V 2 = 0 and in that case, as we will show, the eddy 
feedback f r again falls, but as /3 _1 or for some special non-zonal perturbations 
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(a) 




Figure 3.7: (a) Contours of J r (^,n) in a (0,n) polar plot (n radial and 0 azimuthal) for 

isotropic forcing and ft — 200. This panel shows the vorticity fluxes induced by waves with 
phase lines oriented at an angle 0 in the presence of a non-zonal perturbation with mean flow 
wavenumber n and ip — 15°. Solid (dashed) lines indicate contours with positive (negative) 
values, the contour interval is 2.5 x 10 -3 and the thick line is the zero contour, (b) Locus of 
the roots of T>2(0,n) on the (0, n) plane for non-zonal perturbations with <p = 15°. The roots 
correspond to resonant interaction between waves with phase lines oriented at an angle 0 and 
non-zonal perturbations with mean flow wavenumber n. Thick solid (dashed) lines indicate 
whether the vorticity fluxes produced by the resonant waves are upgradient (downgradient). The 
radial grid interval in both panels is An = 0.25. 


even as /3 -1 / 2 . 

Consider first the emergence of non-zonal structures in the limit /3 1. The 

contribution of each Fourier component of the forcing to the vorticity flux feedback 
T for the case of non-zonal structures at /3 — 200 is shown in Fig. 3.7a. In contrast 
to the cases with /3 <C 1 (or (3 = 0(1), discussed in section 3.4.3), there is only a 
small band of Fourier components that contribute significantly to the vorticity 
flux feedback, as indicated with the narrow tongues in Fig. 3.7a. The reason for 
this selectivity in the response is that for (3 1 the components that produce 

appreciable fluxes, as seen from (3.14), are concentrated on the (0,n) curves that 
satisfy V 2 = 0 (shown in Fig. 3.7b) or equivalently for the (0, n) that satisfy the 
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Figure 3.8: (a) The curves separating the regions in the (n, ip) plane for which V 2 has no roots 
(region D), 2 roots (region B) and four roots (regions A and C). Waves with 0 corresponding 
to two out of the four roots of V 2 found in region A produce upgradient fluxes, (b)-(d) The 
vorticity fluxes T as a function of the angle 6 subtended by the phase lines of the waves and 
the y axis in the presence of a non-zonal perturbation with p = 15° at = 200. The mean flow 
wavenumber is (b) n m 0.25 (in region D), (c) n = 0.5 (in region B), (d) n = 0.592 (in region A) 
and (e) n = 0.75 (in region A). The resonant angles (i.e. the roots of V 2 ) are marked by upper 
(lower) triangles when the waves induce upgradient (downgradient) fluxes. Note that the scale 
in (b) is much smaller. 


resonant condition + uj n = c^k+n (cf. (E.20)). This is the resonant condition 
satisfied when a Rossby wave with wavevector k and frequency forms a resonant 
triad with the non-zonal structure with wavevector n and frequency cj u . We 
concentrate our analysis to these “resonant contributions” because they dominate 
the eddy feedback of non-zonal perturbations for (3 1. 

Resonant triads do not occur for all mean flow perturbations n. For (n, ip) 
in region D of Fig. 3.8a, V 2 has no roots and therefore there are no Fourier 
components with k = (cos 0, sin 6) that form a resonant triad with the mean 
flow perturbation n and the eddy feedback is determined by the sum over the 
non-resonant contributions as illustrated in Fig. 3.8b. In region B of Fig. 3.8a, 
there are only two resonant angles 6. The resonant and non-resonant contribution 
for a typical case in region B is shown in Fig. 3.8c. Note that it is the resonant 
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contributions that determine the eddy feedback. However, they produce a negative 
eddy feedback (a downgradient tendency), which is stabilizing, a result that holds 
for all (n,cp) in region B. In regions A and C, there exist four resonant angles 9 
which dominate the vorticity flux. In C all resonant contributions are stabilizing 
and therefore C is also a stable region. In region A, which at most extends to 
(p = 60° (cf. Appendix F), two of the four resonances give positive contributions 
to f r (cf. Figs. 3.8d,e). Therefore only for (n,(/?) in region A, does a destabilizing 
eddy feedback occur. The largest destabilizing feedback occurs when the positively 
contributing resonances are near coalescence (i.e. as in Fig. 3.8d), which occurs for 
(n, cp) close to the curve separating regions A and B. The reason is that when the 
resonances are apart, as in Figs. 3.8c,e, the significant contributions come from 
near-resonant waves with angles within a band of 0(1//3) around the resonant 
angles and the integrated resonant contributions to the vorticity flux are 0(l//3). 
However, when the resonances are near coalescence, as for the case shown in 
Fig. 3.8d, the band of near-resonant waves contributing significantly increases 
as the integrand assumes a double humped shape and, as shown in Appendix F, 
the destabilizing vorticity flux feedback becomes 0(l/^fj3). Note that as (3 oo, 
the width over which we have significant contributions diminishes and therefore 
f r —y 0 unless an infinite amount of energy is injected exactly at the resonant 
angles (as is assumed in modulational instability studies). 

It can be shown (cf. Appendix F) that the resonant contribution for (3 1 

asymptotically approaches 


N r 


/i R, = ^E 




|p,-I 1 / 2 


(3.17) 


where the subscript j indicates the value of the functions at the j-th out of the 
N r roots of V 2 and p — OqqVz- The values A/}, A),j 5 Pj are all 0(1), whereas rjj is 
always positive and the only quantity that has dependence on /?. It is 0(1) only 
for (n, ip) just above the separating boundaries of regions A and B and regions 
B and D in Fig. 3.8a yielding 1/yfp and is 0( 1/yfff) elsewhere yielding 

fr ~ 1//3, as also qualitatively described above. The sign of the j- th resonant 
contribution to the total eddy feedback depends only on the sign of Mj. For 
(n, ip) just above the boundary separating regions B and D, J\fj < 0 and f r attains 
its minimum value, which corresponds to the largest stabilizing tendency. This 
is illustrated in Fig. 3.9, showing the eddy feedback f r as a function of n. For 
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Figure 3.9: Eddy feedback f r as a function of n for /3 = 200. Positive (negative) values 
correspond to upgradient (downgradient) fluxes. Shown is f r for p = 0° (multiplied by /3 2 ) and 
for p — 15° (multiplied by /3). Also the asymptotic expressions (F.18) for p = 0° and (3.17) for 
p = 15° are shown (dash-dot). The crosses mark the mean flow wavenumbers n = 0.43 and 
n = 0.59 that separate regions A, B and D in Fig. 3.8a for p = 15°. 

(n, y?) just above the boundary separating regions A and B, coalescence of the 
two positive contributing resonances occurs and f r attains its maximum value, 
which corresponds to the largest destabilizing tendency. For small mean flow 
wavenumbers n (corresponding to region D) the eddy feedback is negative and 
0(/3 ~ 2 ) due to the absence of resonant contributions. 

An interesting exception to the results discussed above occurs for the important 
case of zonal jet perturbations (ip — 0 °). In that case, J\fj — 0 in (3.17) as the 
roots of V 2 and J\f coincide and the resonant contribution (3.17) is exactly zero. 
As shown in Fig. 3.10, positive vorticity flux feedback is obtained from a broad 
band of the non-resonant Fourier components with 7 = 0 ~ 0°, corresponding 
to waves with lines of constant phase nearly aligned with the y axis (remember 
that for smaller (3 the region that produces destabilizing fluxes extends up to 
\6\ « 30°). For large /? the vorticity flux f r is always destabilizing for all zonal 
jet perturbations with n < 1, as shown by (F.18) and Fig. 3.9, and the largest 
destabilizing vorticity flux, / r?m ax = (2 + /i)/3 -2 , is obtained for jets with the 
largest allowed scale. The reason for the weak fluxes and the preference for the 
emergence of jets of the largest scale in this limit is understood by noting that 
the stochastically forced eddies for (3 1 propagate with 0(/3) group velocities. 

Therefore in contrast to the limit of f3 1 in which they evolve according to their 
local shear, the forced waves respond to the integrated shear of the sinusoidal 
perturbation over their large propagation extend, which is very weak. 

To summarize: Although zonal jets and most non-zonal perturbations induce 
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Figure 3.10: Contours of .F(0,n) in a (0,n) polar plot (n radial and 0 azimuthal) for zonal 
jet perturbations (ip = 0°) and — 100. Solid (dashed) lines indicate contours with positive 
(negative) values, the contour interval is 2 x 10 —4 , the thick line is the zero contour and the 
radial grid interval is An = 0.25. 

fluxes that decay as 1 //3 2 for large /?, resonant and near resonant interactions 
arrest the decay rate of certain non-zonal perturbations by a factor of 0(/3 3 / 2 ) 
leading to fluxes that decay as 1/ yf]3. This makes the non-zonal perturbations to 
be the most S3T unstable perturbations for /? 1 . Also in contrast to (3 <C 1 

when f r is positive for all n and p (cf. Fig. 3.6), the vorticity flux feedback is 
negative for (n, p) in regions B and D of Fig. 3.8a. As a result, the mean flows 
that produce negative fluxes and are by necessity S3T stable are interestingly in 
the interior of the dumbbell shown in Fig. 3.11, illustrating f r in a polar (n,p) 
plot. The largest destabilizing fluxes occur in the narrow region adjacent to 
the outer boundaries of the dumbbell shape, which demarcates the boundary 
separating regions A and B. Because of the selectivity of the resonances these 
results do not depend on the forcing anisotropy, as we will see in the next section. 

3.4.3 Induced vorticity fluxes for /? ~ 0(1) 

We have seen that in the singular case of isotropic forcing the only process available 
for the emergence of mean flows is the fourth-order anti-diffusive vorticity feedback 
induced by the variation of the group velocity of the forced eddies due to the mean 
flow shear. For /? <C 1, the waves interact with the local shear producing fluxes 
proportional to /? 2 d 4 SU/dy 4 . As f3 increases this growth is reduced since the waves 
interact with an effective integral shear within their propagation extent which is 
weak and eventually, as we have seen in the previous section, for /? 1 the fluxes 

decay as /3~ 2 . Therefore, the fluxes attain their maximum at an intermediate 
value of /?. This occurs for (3 ~ 3.5, as can be seen in Fig. 3.12a where the 
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Figure 3.11: Contours of the eddy feedback f r in a (<p, n) polar plot (n radial and cp azimuthal) 
for the case /3 = 200. Shown are contours of positive values, so the white area corresponds to 
negative values indicating downgradient vorticity fluxes. The contour interval is 10 -3 and the 
radial grid interval is An m 0.25. Note that the feedback factor is always negative (downgradient 
fluxes) for ip > 60° (cf. Appendix F). 
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Figure 3.12: The maximum value of f r over all wavenumbers n for zonal jets (solid), and 
the maximum value of f r over all wavenumbers n and angles ip ^ 0° for non-zonal perturbations 
(dashed) as a function of the planetary vorticity, /3 for the three forcing covariance spectra seen 
in Fig. 3.1 and for p — 1/4. Also shown are the asymptotic expressions (F.4), (F.6) and (F.19) 
(dash-dot) and the slope (dotted). For (i — — 1 zonal jet perturbations are stable for 

ft < 1.67. (b) The mean flow wavenumber n and (c) the angle ip for which the maximum value 
of f r (shown in (a)) is attained. The asymptotes n = l/y/2 (for ft <C 1) and n — 0.5 (for ft 1) 
are shown in (b) (dash-dot) as well as the asymptote ip — 10° (for ft 1) is also shown in (c) 
(dash-dot). 
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Figure 3.13: Contours of the J 7 ^, n) in a (0, n) polar plot (n radial and 0 azimuthal). Shown 
is (a) T for a zonal jet perturbation (<p — 0°) and (c) a non-zonal perturbation with <p — 15° 
when [3 — 2. Panels (b) and (d) are the same as (a) and (c) for the case /3 = 12. In all panels, 
solid (dashed) lines indicate contours with positive (negative) values, the contour interval is 
2 x 10 —3 , the thick lines indicate the zero contour and the radial grid interval is An = 0.25. 
White dashed lines in (c), (d) correspond to the locus of the roots of V 2 (0, n) on the (0, n) plane. 


maximum f r over all (n, cp) is shown. It is demonstrated in the next section that 
this intermediate /? maximizes the S3T instability for all forcing spectra. 

While the eddy-mean flow interaction of both zonal and non-zonal perturbations 
is dominated by the same dynamics when (3 <C 1, for f3 1 the eddy-non-zonal 
flow interaction is dominated by resonances which do not occur for zonal flow 
perturbations. The resonant interactions lead to the possibility of arrested decay 
of the eddy feedback at the rates of /3 -1 / 2 and /3 _1 , instead of the /3~ 2 decay in the 
absence of resonances. The vorticity flux attains its maximum at an intermediate 
value /? ~ 0(1) for non-zonal mean flows as well, which is nonetheless large enough 
for the resonant contributions to reinforce the contribution from the shearing 
mechanism. Figure 3.13 shows the contribution to the eddy feedback induced 
by the various wave components that are excited for two values of /3 (/3 = 2 and 
/? = 12) in the case of zonal jets (cp = 0°) and non-zonal perturbations (cp = 15°). 
As /? increases, the resonant contributions start playing an important role for 
non-zonal perturbations as there is enhanced contribution to the eddy feedback 
in the vicinity of the T >2 — 0 curves, indicated by the white dashed lines. These 
resonant contributions enhance the vorticity fluxes relative to the fluxes obtained 


59 



















for zonal jets and render the non-zonal structures more unstable compared to 
zonal jets when /3 > 3.5 (Bakas and Ioannou, 2014). 

3.5 Effect of anisotropic forcing on S3T instability 

In this section we investigate the effect of the anisotropy of the excitation on 
the S3T instability. The maximum vorticity flux feedback f r for three cases of 
anisotropy (fi = ±1 and fi — 1/4) and for isotropic forcing (/i = 0) is shown in 
Fig. 3.12a. For /3 1, the main contribution to f r for zonal jet perturbations, 

comes from forced waves with nearly meridional constant phase lines (angles near 
0 = 7 = 0°, cf. Fig. 3.10). Therefore, the eddy feedback / r , attains larger (smaller) 
values for a stochastic forcing that injects more (less) power in waves with angles 
near 7 = 0°, that is for positive (negative) anisotropicity factor p (cf. Fig. 3.1). 
The maximum value of f r over all wavenumbers n depends in this case linearly 
on p (cf. Appendix F), 


fr, max — (2 + p) f3 2 + (9(/3 4 ) . (3.18) 

For non-zonal perturbations, the main contribution comes from forced waves 
satisfying the resonant condition + <^n = ^k+n and f r depends only on the 
sum of the resonant contributions. The sign of Afj that determines whether the 
resonant contribution is positive or negative (cf. (3.17)), depends only on the 
sign of sin0j + n/2 and not on the anisotropicity factor (i (cf. (E.19c)). The 
anisotropicity affects only the magnitude of Afj. For any 0 < p < 90° it is found 
that the resonances giving positive contribution occur at angles 6j for which 
IT?I = \0j ~ vl < 45°. A stochastic excitation, which injects more power near 
7 = 0 ° (/i > 0 ) gives larger positive resonant contributions and therefore f r 
increases with fi. However, the effect on the maximum vorticity feedback is weak, 
as the spectral selectivity of the resonances renders the characteristics of the most 
unstable non-zonal structure independent of the spectrum of the forcing. That is, 
the (n, p) that correspond to the maximum f r asymptotes to n ^ 0.5, p « 10° 
(marked with star in Fig. 3.8a) as /? -A oo, a result that is very weakly dependent 
on (i (cf. Figs. 3.12b,c). 

For /? 1, the characteristics of the S3T instability are dependent on the 

anisotropy of the stochastic forcing. The eddy feedback is at leading order 
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proportional to y: 

f r — ^/in 2 ^1 — n 2 ^ cos(2 p) + 0(/3 2 ) . (3.19) 

This shows that there can be upgradient vorticity fluxes leading to S3T instability 
for /? = 0 as long as y cos (2p) > 0. For /i > 0, the maximum f r = /i/32 is 
achieved by zonal jets (p = 0 °), while for fi < 0 any non-zonal perturbation with 
p > 45° can grow, with the maximum f r = \p\/32 achieved for p = 90° when the 
non-zonal perturbations assume the form of jets in the y direction (meridional 
jets) (cf. Fig. 3.12c). 

It is worth noting that Srinivasan and Young (2014) also find that that the 
eddy momentum fluxes are proportional to y when a constant shear flow is 
stochastically forced with power spectrum (3.11). This result is intriguing as the 
two studies address two different physical regimes. This chapter treats the limit 
appropriate for emerging structures in which the shear time is far larger than the 
dissipation time-scale with the fluxes determined by the instantaneous response 
of the eddies on the shear. Srinivasan and Young (2014) study the opposite limit 
in which the mean flow shear is finite and the shear time is much shorter than 
the dissipation time-scale with the fluxes determined by the integrated influence 
of the shear on the eddies over their whole life cycle, which may include complex 
effects such as reflection and absorption at critical levels. 

In summary: 

a. The S3T instability of the homogeneous state is a monotonically increasing 
function of y for all /? 

b. The forced waves that contribute most to the instability are structures 
with small 7 , i.e., waves with phase lines nearly aligned with the y axis, as 
Fig. 3.1a. 

c. The anisotropy of the excitation affects prominently the S3T stability of 
the homogeneous state only for (3 < 3.5. 

3.6 Bibliographical note 

This chapter is an excerpt from the paper by Bakas et al. (2015). The S3T 
instability of the homogeneous turbulent equilibrium was studied by Farrell and 
Ioannou (2007). Analytical results for zonal jet perturbations for /? = 0 and finite 
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doubly periodic domains were obtained by Bakas and Ioannou (2011). Results 
for infinite /3-planes and for any /3 were obtained by Srinivasan and Young (2012). 
The forcing spectrum used in this chapter was introduced by Srinivasan and 
Young (2014). The dispersion relation for the stability of non-zonal perturbations 
was derived by Bakas and Ioannou (2013a). A physical interpretation of the S3T 
instability of the homogeneous equilibrium to zonal jet perturbations for small /3 
and n is discussed in Bakas and Ioannou (2013b). 
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A note on non-dimensional units used in the 

following chapters 


In the next chapters we will not scale our fields and parameters as described 
in (2.3). Instead, we will non-dimensionalize everything using typical values 
that correspond to the Earth’s midlatitude atmosphere, that is a length scale of 
L = 5000 km and a velocity of U — 40 ms -1 . Using this scales the time unit 
is T — 1.5 day and the Earth’s meridional planetary vorticity gradient at the 
midlatitudes corresponds to the non-dimensional value /? = 10. 

All numerical simulations that will be presented in the following chapters will be 
implement using periodic boundary conditions on a /3-plane with non-dimensional 
size L x — L y — 2n. Therefore non-dimensional wavenumbers assume only integer 
values. 
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4 

Relation of the S3T system with the 4MT 
system of the modulational instability of 

Rossby waves 


In this chapter we investigate the relation of the modulational instability (MI) of 
Rossby waves with the S3T theory regarding the emergence and equilibration of 
large-scale structures in /3-plane turbulence. It was established by Lorenz (1972) 
and Gill (1974) that Rossby waves are hydrodynamically unstable, and under 
certain conditions the greatest instability is a zonal jet. This instability is an 
instability that has been characterized in the literature as a MI because of its 
similarity with the Benjamin-Feir instability of surface gravity waves (Benjamin, 
1967; Yuen and Lake, 1980). More recently this instability has been proposed to 
be the mechanism for the formation of zonal jets in barotropic but also baroclinic 
turbulence (Berloff et al., 2009; Connaughton et al., 2010), in the sense that at 
the Rhines’s scale the turbulent state is dominated by relatively coherent wave 
structures that become modulationally unstable and give rise to jets. In this 
chapter we demonstrate the formal equivalence between the 4MT system, that 
approximates well the MI of coherent Rossby waves, and the S3T instability of a 
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homogeneous turbulent state that has the power spectrum of the Rossby wave 
that undergoes ML This equivalence embeds the MI results into the more general 
and physical framework of S3T which can address the instability of more general 
states, like the structural instability of the attractor of a turbulent flow. We also 
compare the predictions of the 4MT and S3T systems with nonlinear simulations 
regarding the initiation of the MI and its equilibration. We demonstrate that 
the 4MT dynamical framework is inadequate for capturing the finite amplitude 
equilibration of the instabilities. 

4.1 MI of a Rossby wave and the 4MT approximation 

Consider the stability of a Rossby wave with streamfunction ^ p = A cos(p-x— u p t) 
(and vorticity £ p = —p‘ 2r tp p ) that satisfies the inviscid barotropic vorticity equation: 

dtCp + J (V’P) Cp + P ■ x) = o . (4.1) 

The stability of these nonlinear traveling wave solutions, referred to in MI studies 
as the primary waves, is addressed by perturbing the primary wave , i.e., by 
writing £ = £ p + S( and studying the evolution of the perturbation S( in the 
linear approximation, 

d t 6( = £(( p )6( , (4.2) 

where £(£ p ) is a time-dependent linear operator. With the change of the frame 
of reference: 

(z x (3) t 

X 0 = X-«- 

p z 

the primary wave assumes the stationary form: ^ p = A cos(p • xq), and the 
operator C becomes time-independent but with the spatial periodicity of the 
primary wave. The eigensolutions of (4.2) according to Bloch’s theorem, are 

^Cn(x 0 , t) = e Snt e in ' x ° g(x 0 ) , (4.4) 

and each eigenfunction is indexed by a wavevector n which satisfies |n| < |p|/2. 
The function g is a periodic function with the periodicity of the primary wave i/j p 
(cf. Appendix D, eq. (D.7)), and can be assumed in the form: 

+oo 

<7(xo)= E «me imp ' x °. (4.5) 
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In the original coordinates, because: 


+oo 

g(* o) -»• g (x - (z x /3) t/p 2 ) = X) e im (p- x -p 4 ) , 

m =—oo 

gin-x 0 gin- [x—(zx/3) t/p 2 ] _ e i(l-n 2 /p 2 )w n £ e i(n-x-u; n t) 


(4.6a) 

(4.6b) 


(c^k denotes the frequency of a Rossby wave k, cf. (3.6)), the eigenfunction n can 


be written as 

<Kn(x, t) = e Sny(l-n 2 / P >n t^ n (x, t) , (4.7) 


with 


+oo 

5(n(x,t) = a 0 e i(n ' x_aJnt) + 51 e i[(n+mp)'X-k+m,p)i] _ (4. 8 ) 

m =—oo 
m/0 

Written in this form the eigenfunction (4.8) is a superposition of a nonlinear Rossby 
wave solution of (4.1) (the ao harmonic) and satellite modes with wavenumbers 
n zb mp, m — 1,2,..., that are Rossby wave solutions only when cj n + mu p = 

^n+mp- 

By inserting (4.4)-(4.5) into (4.2) an infinite homogeneous linear system for the 
coefficients a m is obtained. The eigenvalues s n are obtained from the requirement 
that this system has non-trivial solutions. This implies that the s n are roots of 
the associated characteristic polynomial, which is nominally of infinite degree. 
However, because the physically realizable solutions correspond to the convergent 
series (4.8), the coefficients of physically realizable eigenfunctions will have the 
property that a m 4 0 as m 4 zboo; actually a m ~ b m /m\ for some constant 
b (Lorenz, 1972). This enables us to determine accurate approximations of the 
eigenvalues from finite truncations of this infinite system. One obtains a good 
approximation of the eigenvalue even if only terms up to \m\ = 1 are kept. This 
truncation is referred to as the 4 mode truncation, or “4MT” system, because 
only four waves are allowed to interact: the primary wave p and perturbation 
waves n, nip. It can be also shown from a Fjprtoft type argument that unstable 
eigenvalues with s r = Re(s) > 0 exist only for |n| < |p| (Lorenz, 1972). The 
instability manifests as a modulation of the amplitude, as shown in Fig. 4.1. 
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Figure 4.1: Modulational instability of a primary wave f(x, t) = cos (kx — Ukt ) with k = 10 
to a large-scale perturbation Sf(x, t) that grows at rate s r ' Sf(x, t) = 0.1 e Srt cos (nx — u n t) with 
n — k/ 5. Solid lines correspond to / + Sf while dash-dotted lines correspond to the unperturbed 
/. Initially the “energy” of large-scale perturbation, E m , is 1% of the “energy” of the finite 
amplitude of the primary wave, E p . The instability manifests as a modulation to the wave 
amplitude of the primary wave (cf. panel (b)). 


4.2 Equivalence of the MI in the 4MT approximation and the S3T 

STABILITY OF A HOMOGENEOUS TURBULENT STATE 

There is a close relation between the 4MT approximation of the MI and the 
S3T. Parker and Krommes (2015) have shown that in the inviscid limit there is 
a formal equivalence between the modulational instability of the Rossby wave, 
= A cos(p-x — Upt), in the 4MT approximation with the S3T instability of the 
homogeneous state with eddy vorticity covariance with the same power spectrum 
as the Rossby wave, i.e., with (7 e (k) = (27r) 2 p 4 |A| 2 [5(k — p) + S( k + p)]. The 
connection is formal because physically the two problems are very different. In 
MI the stability of a basic state in the form of a coherent Rossby plane wave is 
studied, while S3T addresses the statistical stability of an incoherent state with 
equilibrium covariance having the power spectrum of the Rossby wave. In that 
sense, as noted by Parker and Krommes (2015), S3T stability analysis embeds 
the modulational instability results into a more general physical framework. 

We proceed here to show that this result does not only hold for monochromatic 
waves but can be generalized to any solution of the barotropic vorticity equation. 
That is, we show the formal equivalence between the MI of any time-dependent 
solution of the barotropic equations with stationary power spectrum in the dynam¬ 
ical framework of a generalized 4MT with the S3T instability of a homogeneous 
state with the same power spectrum. The proof can be found in Appendix G. 
Such a nonlinear solution of the inviscid barotropic vorticity equations is for 
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example a superposition of any number of Rossby waves: 

N 

ip = A J cos (Pj • x “ u Pj t) , (4.9) 

j =1 v ----—' 

Ip j\=p ^ p j 

all with the same total wavenumber 1 , that forms a non-dispersive structure moving 
westwards (cf. Appendix G). The MI of (4.9) can be carried out in a similar 
manner as described in section 4.1. The generalized 4MT system is obtained 
by keeping the N primary waves, ip Pj , the perturbation wave, e in x °, and the 
corresponding 2N satellite modes, ^)' x °. The eigenvalue relation in this 

truncation coincides with the S3T eigenvalue relation of the equilibrium covariance 
with spectral power: 


C e (k) 


N 

( 2 vr ) V ^|^| 2 


3 = 1 


5(k — p j) + 5( k + p j) 


(4.10) 


The MI of a base state in the form of (4.9) as well as the mechanisms responsible 
for instability have been studied by Lee and Smith (2003). 


4.3 Comparison of MI and S3T predictions with nonlinear simula¬ 
tions FOR THE EMERGENCE AND EQUILIBRATION OF JETS 

Connaughton et al. (2010) compared the predictions of the 4MT system with 
direct numerical simulations and found that the 4MT system captures the initial 
growth of the instability, but fails to predict the later stages of zonal flow evolution. 
Contrary to the 4MT system, S3T dynamics capture both the emergence of the 
large-scale flow instability and also its equilibration. Here we present an example 
of jet emergence and equilibration as predicted by the 4MT and S3T systems and 
compare them with nonlinear simulations of the barotropic vorticity equation. 
Details regarding the methods used for performing the numerical simulations can 
be found in Appendix C. 

We start by performing a simulation of the inviscid and unforced version of 
the NL system (2.1), referred to as NLi nv . We initiate the simulation with a state 
^(x,£ = 0) that consists of a primary wave — Acos(p • x) with p = (7,0) 
and energy E p (t = 0) = 2 x 10 -3 and a zonal jet perturbation ip n = a cos (n • x) 

1 The vorticity of (4.9), Aip = —p 2 ^, is proportional to ip and as a result J{ip, Aip) = 0. 
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with n = (0,3) and energy E m {t = 0) = 10 -9 — 0.5 x 10 -7 i^(£ = 0). This 
n = (0,3) zonal flow perturbation has been chosen because it is predicted by 
both S3T or 4MT to be the most unstable large-scale structure. We also perform 
a 4MT simulation which is initiated with the initial state of the NLi nv . In the 
4MT dynamics we allow only interactions between the Fourier components with 
wavenumbers =tp, ±n and =L(nd=p). The evolution of the zonal mean flow energy, 
E m , in the two simulations is plotted in Fig. 4.2; for comparison we also plot the 
growth of the energy of the emerging instability as predicted by S3T (or 4MT). 
Snapshots of the evolution of the flow streamfunction, are shown in Fig. 4.3 
(NLi nv ) and Fig. 4.5 (4MT), for the time instants marked in Fig. 4.2. Additionally, 
we perform an integration of the stochastically forced-dissipative NL system (2.1), 
in which the forcing excites structures cos (p • x + 9) with 6 a randomly chosen 
phase. (The spatial covariance of this forcing is Q(x a — x&) ~ cos [p • (x 0 — x&)].) 
For the chosen coefficient of linear damping, r, the energy input rate, e, is 
adjusted so that the steady state equilibrium total energy of the stochastically 
forced flow is equal to the total energy of the primary wave of NLi nv and 4MT, 
i.e., ^ = 2 rE p (t = 0) (cf. (A.22)). The evolution of the zonal mean flow as 
well as snapshots of the flow streamfunction are shown in Fig. 4.2 and Fig. 4.4 
respectively. 
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Figure 4.2: Evolution of the zonal energy E rn (/c cc = 0) in the 4MT system (solid red line), 
the NLinv system (dash-dotted green line) and the forced-dissipative NL system (solid black 
line). Both the NLi nv and 4MT systems are initiated with a plane wave with wavenumber 
p = (7, 0) and E p (t = 0) = 2 x 10 -3 and a zonal jet perturbation with wavenumber n = (0, 3) 
and energy Emit — 0) = 10 -9 . The parameters for the NL are: linear damping coefficient 
r — 0.01, stochastic forcing with single harmonics with wavenumber p and energy injection rate: 
s = 2rE p (t = 0) = 4 x 10 -5 . The predicted growth of the n = (0, 3) zonal jet perturbation by 
S3T is shown with the dashed line. Remarkably, the energy E m (t) of the mean flow grows at the 
same rate in the unforced and inviscid NLi nv and the forced-dissipative NL. Typical snapshots 
of the streamfunction fields for the three simulations are shown in Fig. 4.3 (NLi nv ), Fig. 4.4 
(NL) and Fig. 4.5 (4MT), for the times marked with circles. Both the S3T and 4MT predict the 
initial growth of the mean flow but the 4MT fails to capture the finite amplitude state of the 
system. In all simulations = 4.9. 
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(a) ip(*,t = 0) (b) r(x. t -= 18) ( c ) V’(x,i = 64) 



Figure 4.3: Snapshots of streamfunction -0(x, t) together with the zonally averaged zonal 
velocity U ( y , £) (thick black line) for the NLi nv system at the indicated times with circles in 
Fig. 4.2. Initially the zonal mean perturbation n = (0,3) grows to finite amplitude (panels 
(a)-(c)) and at t ~ 200 the zonal flow reorganizes and becomes a (1,4) westward traveling wave 
(panels (d)-(f)). 
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Figure 4.4: Same as Fig. 4.3 but for the stochastically forced-dissipative NL system. 
Remarkably, the NL system exhibits the same large-scale structure evolution with the NLi nv 
and transitions at approximately the same time to a (1,4) traveling wave structure. 
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Figure 4.5: Snapshots of streamfunction ^(x, t) together with the zonally averaged zonal 
velocity U(y,t) (thick black line) for the 4MT system at the indicated times with circles in 
Fig. 4.2. Initially the zonal mean perturbation n s= (0,3) grows to finite amplitude but then it 
alternates between a state with strong zonal mean flow component (i.e. panel (d)) and a state 
with weak zonal flow component (i.e. panel (e)). 
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Initially the zonal mean flow E m (t ) in the NLi nv and NL grows at the rate 
predicted by the S3T and 4MT stability analysis. In both the NL and NLi nv the 
amplitude of the zonal flow reaches a plateau and then the flow reorganizes at 
t ~ 220 producing a traveling wave moving westwards with maximal power at 
wavenumber (1,4) (cf. Fig. 4.3e,f). Remarkably, the forced-dissipative NL under¬ 
goes the same flow reorganization to a traveling wave mean flow at approximately 
the same time as the NLi nv (cf. Figs 4.2 and 4.4). The 4MT system however, fails 
to capture this structural reorganization of the flow. Instead, it oscillates between 
a state with a strong zonal mean flow component (i.e. Fig. 4.5d) and a state with 
weak zonal mean flow component (i.e. Fig. 4.5e) that have no reflection in the 
NLi nv . 

In order to investigate whether the S3T system is able to produce the NL large- 
scale flow state we perform a forced-dissipative S3T time integration of (2.13) 
with the parameters of the NL simulation. Snapshots of the large-scale flow 
streamfunction, T, that emerges in the S3T simulation are plotted in Fig. 4.7. 
First a zonal mean flow emerges, the zonal flow equilibrates producing finite 
amplitude jets, which then become unstable and give way to a traveling wave 
with structure similar to that of the NL simulation (see Fig. 4.7e,f). 

At this point we want to emphasize that the S3T that succeeded to faithfully 
produce the NL flow state was the S3T system (2.13) in which the ensemble mean 
was identified with an average over fast time scales. The simplest S3T system 
in which the ensemble mean is identified with a zonal mean, that we denote 
as S3Tz and obeys equations (2.21), is able to reproduce the initial instability 
and equilibration of the zonal jet but is incapable to capture the transition to a 
non-zonal large-scale flow and instead it equilibrates to a zonal jet state with 3 
jets (see Fig. 4.8). A comparison of the evolution of the zonal mean flow energy 
for the S3T and S3Tz systems is shown in Fig. 4.6. The S3Tz system mean flow 
energy evolution initially coincides with the S3T energy evolution and at t « 150 
the two energy evolutions diverge: the S3Tz system is attracted to a zonal jet 
mean flow while the S3T system transitions to a turbulent state characterized by 
a traveling wave mean flow with most power at wavenumber (1,4). 

The stationary statistical equilibrium that the S3Tz system is attracted is also 
a stationary state of the S3T system and it can be shown that is unstable to 
non-zonal mean flow perturbations. With the methods described in chapter 6 we 
determine that the most unstable eigenfunction of the large-scale flow corresponds 
to a traveling wave with wavenumber (1,4) (cf. section 6.2). This demonstrates 
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Figure 4.6: Evolution of the zonal energy Em(k x = 0) in logarithmic scale (panel (a)) and in 
linear scale (panel (b)) in the 4MT system (solid red line), the NL system (black red line), the 
S3T system (blue solid line) and the S3Tz system (dashed line). It can be seen that the S3T 
prediction follows closely the NL. The S3Tz system mean flow energy evolution initially coincide 
with the S3T evolution and at t « 150 the two systems diverge. The S3Tz system equilibrates to 
a zonal mean flow statistical equilibrium while the S3T system transitions to a traveling wave 
mean flow with most power at wavenumber (1,4). For all integrations the planetary vorticity 
gradient is (3 = 4.9. Circles mark the time instants for the snapshots plotted in Figs. 4.4, 4.5, 4.7 
and 4.8. 


that the final state of the NL is an equilibrated secondary instability of the S3T 
dynamics. 


74 

















(a) *(x,t = 0) (b) *(x,t = 48) 


6 



6 


4 



4 

) 





\ 

2 



2 

/ 

0 

, 

, 

0 

1) . 


X 


X 


(c) ®(x, t = 64) 

5 


5 


(d) *(x, t = 224) ( e ) tt(x, t* 292) (f) (x, t:* 772) 



6 

4 

!i ^ 

6 

4 



2 

i / 

2 

- 






n 

_1_1_1_L 

0 

i\ 

__1_I_ ^ _L 

0 

, y 


024602460246 

XXX 


0.02 

0 


- 0.02 


Figure 4.7: Snapshots of mean flow streamfunction T(x, t) together with the zonally averaged 
zonal velocity U(y,t) (thick black line) for the S3T system at the indicated times with circles in 
Fig. 4.6b. Initially a zonal mean flow emerges and equilibrates to finite amplitude (panel (d)). 
This state however becomes unstable and transitions to a traveling wave with structure similar 
to that of the NL. 
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Figure 4.8: Snapshots of mean flow streamfunction T(y, t) together with the zonally averaged 
zonal velocity U(y,t) (thick black line) for the S3Tz system at the indicated times with circles 
in Fig. 4.6b. In the S3Tz system the mean flow is by construction x -independent and therefore 
is able to reproduce the initial instability and equilibration of the zonal jet but is incapable to 
capture the transition to a non-zonal large-scale flow. Instead the zonal flow initially grows and 
then equilibrates to a zonal jet state with 3 jets. 


75 
























































5 

Emergence and equilibration of jets in 
/3-plane turbulence as predicted by S3T and 
its reflection in nonlinear simulations 


5.1 Introduction 

Stochastic structural stability theory (S3T) addresses turbulent jet dynamics as a 
two-way interaction between the mean flow and its consistent field of turbulent 
eddies (Farrell and Ioannou, 2003). The mean flow is supported by its interaction 
with a broad turbulence spectrum through non-local interactions in wavenumber 
space. In fact, S3T is a non-equilibrium statistical theory that provides a closure 
comprising a dynamics for the evolution of the mean flow together with its 
consistent field of eddies. In S3T the dynamics of the turbulence statistics 
required by this closure are obtained from a stochastic turbulence model (STM), 
which provides accurate eddy statistics for the atmosphere at large scale (Farrell 
and Ioannou, 1993, 1994, 1995; Zhang and Held, 1999). 

Marston et al. (2008) have shown that the S3T system is obtained by truncating 
the infinite hierarchy of cumulant expansions to second order and they refer to 
the S3T system as the second order cumulant expansion (CE2). In S3T, jets 
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initially arise as a linear instability of the interaction between an infinitesimal 
jet perturbation and the associated eddy field and finite amplitude jets result 
from nonlinear equilibria continuing from these instabilities. Analysis of this jet 
formation instability determines the bifurcation structure of the jet formation 
process as a function of parameters. In addition to jet formation bifurcations, S3T 
predicts jet breakdown bifurcations as well as the structure of the emergent jets, 
the structure of the finite amplitude equilibrium jets they continue to, and the 
structure of the turbulence accompanying the jets. Moreover, S3T is a dynamics 
so it predicts the time dependent trajectory of the statistical mean turbulent 
state as it evolves and, remarkably, the mean turbulent state is often predicted 
by S3T to be time dependent in the sense that the statistical mean state of the 
turbulence evolves in a manner predicted by the theory (Farrell and Ioannou, 
2009b). The formation of zonal jets in planetary turbulence was studied as a 
bifurcation problem in S3T by Farrell and Ioannou (2003, 2007, 2008, 2009a, c); 
Bakas and Ioannou (2011); Srinivasan and Young (2012); Parker and Krommes 
(2014a). A continuous formulation of S3T developed by Srinivasan and Young 
(2012) has facilitated analysis of the physical processes that give rise to the S3T 
instability and construction of analytic expressions for the growth rates of the 
S3T instability in homogeneous /3-plane turbulence (Srinivasan and Young, 2012; 
Bakas and Ioannou, 2013b; Bakas et al., 2015). Recently, the analogy between 
the dynamics of pattern formation and zonal jet emergence in the context of S3T 
was studied by Parker and Krommes (2014a). 

Relating S3T to jet dynamics in fully nonlinear turbulence is facilitated by 
studying the quasi-linear (QL) model which is intermediate between the nonlinear 
model and S3T. The QL approximation to the full nonlinear dynamics (NL) results 
when eddy-eddy interactions are not explicitly included in the dynamics but are 
either neglected entirely or replaced with a simple stochastic parameterization, so 
that no turbulent cascade occurs in the equations for the eddies, while interaction 
between the eddies and the zonal mean flow is retained fully in the zonal mean 
equation. S3T is essentially QL with the additional assumption of an infinite 
ensemble of eddies replacing the single realization evolved under QL. Although the 
dynamics of S3T and QL are essentially the same, by making the approximation 
of an infinite ensemble of eddies, the S3T equations provide an autonomous and 
fluctuation-free dynamics of the statistical mean turbulent state, which transforms 
QL from a simulation of turbulence into a predictive theory of turbulence. 

A fundamental attribute of QL/S3T is that the nonlinear eddy-eddy cascade 
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of NL is suppressed in these systems. It follows that agreement in predictions 
of jet formation and equilibration between NL and QL/S3T provides compelling 
evidence that cascades are not required for jet formation and theoretical support 
for observations showing that the turbulent transfers of momentum maintaining 
finite amplitude jets are non-local in spectral space. 

Previous studies demonstrated that unstable jets maintained by mean flow 
body forcing can be equilibrated using QL dynamics (Schoeberl and Lindzen, 
1984; DelSole and Farrell, 1996; O’Gorman and Schneider, 2007; Marston et ah, 
2008). In contrast to these studies, in this work we investigate the spontaneous 
emergence and equilibration of jets from homogeneous turbulence in the absence 
of any coherent external forcing at the jet scale. S3T predicts that infinitesimal 
perturbations with zonal jet form organize homogeneous turbulence to produce 
systematic up-gradient fluxes giving rise to exponential jet growth and eventually 
to the establishment of finite amplitude equilibrium jets. Specifically, the S3T 
equations predict initial formation of jets by the most unstable eigenmode of 
the linearized S3T dynamics. In agreement with S3T, Srinivasan and Young 
(2012) found that their NL simulations exhibit jet emergence from a homogeneous 
turbulent state with subsequent establishment of finite amplitude jets, while 
noting quantitative differences between bifurcation parameter values predicted 
by S3T and the parameter values for which jets were observed to emerge in 
NL. Tobias and Marston (2013) also investigated the correspondence of CE2 
simulations of jet formation with corresponding NL simulations and found that 
CE2 reproduces the jet structure, although they noted some differences in the 
second cumulant, and suggested a remedy by inclusion of higher cumulants. 

In this chapter we use NL and its QL counterpart together with S3T to examine 
further the dynamics of emergence and equilibration of jets from turbulence. 
Qualitative agreement in bifurcation behavior among these systems, which is 
obtained for all the spatial turbulence forcing distributions studied, confirms that 
the S3T instability mechanism is responsible for the formation and equilibration 
of jets. Quantitative agreement is obtained for bifurcation parameters between NL 
and QL/S3T when account is taken of the modification of the turbulent spectrum 
that occurs in NL but not in QL/S3T. Remarkably, a primary component of this 
spectral modification can itself be traced to S3T instability, but of non-zonal 
rather than of zonal form. We investigate the formation and equilibration of these 
non-zonal S3T instabilities and the effect these structures have on the equilibrium 
spectrum of /3-plane turbulence. We also investigate circumstances under which 


78 


non-zonal structures are modified and suppressed by the formation of zonal jets. 

A dynamic of potential importance to climate is the possibility of multiple 
equilibria of the statistical mean turbulent state being supported with the same 
system parameters (Farrell and Ioannou, 2003, 2007; Parker and Krommes, 2014a). 
We verify existence of multiple equilibria, predicted by S3T, in our NL simulations. 
Finally, we show that weak jets result from stochastic excitation by the turbulence 
of stable S3T modes, which demonstrates the physical reality of the stable S3T 
modes. Turbulent fluctuation induced excitation of these weak local jets and 
the weak but zonally extended jets that form at slight supercriticality in the 
jet instability bifurcation may explain the enigmatic latent jets of Berloff et al. 
( 2011 ). 

Since the emergence and equilibration of jets is addressed throughout this 
chapter the zonal mean-eddy decomposition for the flow fields (2.14) is used. 
Therefore the NL, QL and S3T dynamics discussed in this chapter use the zonal 
mean-eddy decomposition which was presented in section 2.2. The NL system is 


d t U = v'(' ~rU , 

dtC = MU) C' + dy (/C 7 ) - V - (u'CO +Ve£ , 

with 

MU) = ~Ud x -{p- dl y U) S.A- 1 - r , 

the QL system is 

d t U = vp -rU , 

dtC = M u ) C + , 


and the S3T system (S3Tz) is 

1 


d t U = 


l dx a +\ lu x b )^ab 

X a =x b 

9tC a b = v4z ,a(U) + A z ,b(U) C a b + sQab ■ 


-rU , 


(5.1a) 

(5.1b) 

(5.2) 

(5.3a) 

(5.3b) 

(5.4a) 

(5.4b) 


Throughout this chapter where we mention S3T we refer to the zonal mean- 
eddy decomposition S3Tz system. The S3Tz system does not allow for mean flows 
with non-zonal structure and therefore S3T instabilities within (5.4) are only 
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zonal jet perturbations. However, as it is demonstrated in sections 5.4 and 5.5, 
the emergence of non-zonal coherent structures is of great importance in the 
quantitive predictions of S3Tz for jet emergence and equilibration. 

5.2 Specification of the stochastic forcing structure 

Because the S3T instability mechanism that results in jet bifurcation from a 
homogeneous turbulent state differs for isotropic and non-isotropic turbulence 
(cf. chapter 3), we consider examples of both isotropic and non-isotropic turbulence 
forcing. The jet forming instability in the case of homogeneous, non-isotropic 
forcing arises from the up-gradient fluxes induced by shearing of the turbulence 
by the infinitesimal perturbation jet, while the up-gradient fluxes for the case 
of homogeneous isotropic forcing arise from the refraction of the eddies caused 
by the variation in the potential vorticity gradient induced by the infinitesimal 
perturbation jet. 

Three stochastic forcing structures will be used in our investigation of the 
correspondence among S3T, QL and NL dynamics. The first independently excites 
a set of zonal wavenumbers. This stochastic forcing is spatially homogeneous 
but not isotropic and will be denoted as NIF (non-isotropic forcing). The second 
forcing, denoted IRFn, is an isotropic narrow ring forcing concentrated near a 
single total wavenumber, kf. The third forcing we use, denoted IRFw, is an 
isotropic ring forcing in which the forcing is distributed over a wide annular region 
in wavenumber space around the central total wavenumber. Specification of these 
stochastic forcing structures are given in Appendix H. Plots of the corresponding 
forcing covariance power spectra together with instantaneous realizations both in 
vorticity and streamfunction for the three types of forcing structures are shown in 
Fig. 5.1. Note, that the IRFn ring forcing is peculiar in that it primarily excites 
vortices of scale 1 /kf that are evident in both the vorticity and streamfunction 
fields, while IRFw produces a streamfunction field dominated by large scale 
structure similar to the fields excited by the other broadband forcings. 

5.3 Bifurcations predicted by S3T and their reflection in QL and 
NL SIMULATIONS 

We examine the counterpart in NL and QL simulations of the S3T structural 
instability by comparing the evolution of the domain averaged energy of the zonal 
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Figure 5.1: Contour plots of the spatial Fourier coefficients of the forcing vorticity covariances, 
Qk (cf. (2.12)), used in this study and example realizations of the forcing. Panel (a): Q k for 
NIF with zonal wavenumbers k x = 1,..., 14 and d = 1/5. Panel (d): Q k for IRFn at kf = 14 
and 5kf = 1. Panel (g) Q k for IRFw at kf = 14 and 5k / = 8/V2. In (b), (e) and (h) are shown 
realizations of these forcings in the vorticity field, and in (c), (f) and (i) are shown realizations 
in the streamfunction field. 


flow: 

E m {t) = -d— (\u 2 d 2 x. (5.5) 

-L'x-L'y J z 

The amplitude of the zonal flow is measured with the zonal mean flow index (zmf) 
defined as zmf = + E p ), where E m is the time average of the domain 

averaged zonal mean flow, given in (5.5), and E p is the time average of the domain 
averaged kinetic energy of the eddies, 

E p(t) = y—j— [ \W\ 2 d 2 x . (5.6) 

Zmf is shown as a function of the energy input rate in Fig. 5.2a for NIF 
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Figure 5.2: Bifurcation structure comparison for jet formation in S3T, QL, and NL. Shown 
is the zmf index of jet equilibria for (a) NIF and (b) for IRFn forcing as a function of the forcing 
amplitude £/£ CjZ for the NL simulation (dash-dot and circles), the QL simulation (dashed and 
dots) and the corresponding S3Ta simulation (solid). The bifurcation diagram and the structure 
of the jet agree in the QL and S3Ta simulation, but the bifurcation in the NL simulations occurs 
at ~ lle c ,z for NIF and at at £c^ L ~ 4s CjZ for IRFn. Agreement between NL and S3T 

predictions is obtained if the S3T is forced with the spectrum that reflects the modification 
of the equilibrium NIF or IRFn spectrum respectively by eddy-eddy interactions (the results 
of this S3T simulation is indicated as S3Tb, see discussion at section 5.4). (For IRFn this 
spectrum is shown in Fig. 5.5c.) This figure shows that the structural stability of jets in NL 
simulations is captured by the S3T if account is taken of the nonlinear modification of the 
spectrum. Parameters: /3 = 10, r — 0.01. 


forcing and in Fig. 5.2b for IRFn forcing with r — 0.01. The fundamental 
qualitative prediction of S3T that jets form as a bifurcation in the strength of 
the turbulence forcing is verified in these plots. Agreement in the value of the 
bifurcation parameter is also obtained between S3T and QL while the bifurcation 
parameter is substantially larger in NL. For example, the NL simulations bifurcate 
at ~ ll£ c , z under NIF forcing and at ~ 4s c?z under IRFn forcing. 

Similar behavior was noted by Srinivasan and Young (2012). The reason for 
this difference between the NL and S3T bifurcation curves is revelatory of the 
underlying dynamics of the bifurcation, as we explain in section 5.4. 

S3T dynamics not only predicts the emergence of zonal jets as a bifurcation 
in turbulence forcing, but also predicts the structure of the finite amplitude jets 
that result from equilibration of the initial jet formation instability. These finite 
amplitude jets correspond to fixed points of the S3T dynamics. An example for 
IRFn strongly forced with e — 100s c ^ z and with damping r = 0.01 is shown in 
Fig. 5.3. This example demonstrates the essential similarity among the jets in 
NL, QL and S3T simulations. 

Under strong turbulence forcing the initial S3T jet formation instability typically 
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Figure 5.3: Hovmoller diagrams of jet emergence in NL, QL and S3T simulations with IRFn 
forcing at energy input rate £ = 100£ c ,z- Shown is U(y,t) for the NL (panel (a)), QL (panel 
(b)) and S3T (panel (c)) simulations. Also shown are the equilibrium jets (panel (d)) in the NL 
(dash-dot), QL (dashed), and S3T (solid) simulations. There is very good agreement between 
the jet structure in the NL, QL and S3T simulations, despite the difference in the zmf index 
among them (cf. Fig. 5.2b). Moreover, in all three simulations similar jet mergers are observed, 
leading eventually to final equilibrium jets with smaller meridional wavenumber than that of the 
initial instability. Parameters are [3 — 10, r = 0.01. 


reaches final equilibrium as a finite amplitude jet at a wavenumber smaller than 
that of the initial instability. An example is the case of IRFn at e = 100£ C5Z 
shown in Fig. 5.3. In this example, the jets emerge in S3T initially with zonal 
wavenumber n y = 10, in agreement with the prediction of the S3T instability of 
the homogeneous equilibrium, but eventually equilibrate at wavenumber n y = 3 
following a series of jet mergers, as seen in the Hovmoller diagram. Similar 
dynamics are evident in the NL and QL simulations. This behavior can be 
rationalized by noting that if the wavenumber of the jet remains fixed then as 
jet amplitude continues to increase under strong turbulence forcing violation of 
the Rayleigh-Kuo stability criterion would necessarily occur. By transitioning to 
a lower wavenumber the flow is able to forestall this occurrence of inflectional 
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instability. However, detailed analysis of the S3T stability of the finite amplitude 
equilibria near the point of jet merger reveals that these mergers coincide with the 
inception of a structural instability associated with eddy-mean flow interaction, 
which precedes the occurrence of hydrodynamic instability of the jet (Farrell and 
Ioannou, 2003, 2007). 1 The stability of finite amplitude S3T equilibria will be 
discussed in chapter 6. 



Figure 5.4: (a) Hovmoller diagram showing details of the jet mergers for t < 350 in the 
S3T simulation in Fig. 5.3. In (b) is shown the amplitude of the jet maxima that appear in 

(a) . Note that only the prograde jets merge. The bottom panels show the mean potential 
vorticity gradient p — U yy as a function of y at the times indicated by vertical lines in (a) and 

(b) . These graphs show that the structure of the jets is configured at each instant to satisfy the 
Rayleigh-Kuo stability criterion and that jet mergers are the mechanism in S3T for avoiding 
inflectional instability. Decrease in the amplitude of the jets prior to merger indicates increased 
downgradient vorticity fluxes as the flow approaches hydrodynamic neutrality. 

x Jet mergers occur in the Ginzburg-Landau equations that govern the dynamics of the S3T 
instability of the homogeneous equilibrium state for parameter values for which the system 
is close to marginal stability (Parker and Krommes, 2014a). However, these mergers in the 
Ginzburg-Landau equations are associated with equilibration of the Eckhaus instability, rather 
than equilibration of the inflectional instability associated with violation of the Rayleigh-Kuo 
criterion as is the case for mergers of finite amplitude jets (cf. Fig. 5.4). Characteristic of this 
difference is that in the case of the Ginzburg-Landau equations both the prograde and retrograde 
jets merge, while in the case of the finite amplitude jets only the prograde jets merge. The same 
phenomenology as in the Ginzburg-Landau equations occurs in the case of the Cahn-Hilliard 
equations that govern the dynamics of marginally stable jets in the modulational instability 
study of Manfroi and Young (1999). 
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5.4 Influence of the turbulence spectrum on the S3T jet formation 

INSTABILITY 

Both QL and S3T dynamics exclude interactions among eddies and include 
only the non-local interactions between jets, with k x — 0, and eddies, with 
k x ^ 0. Therefore, there is no enstrophy or energy cascade in wavenumber 
space in either QL or S3T dynamics and the homogeneous S3T equilibrium state 
(cf. (2.22)) has spectrum, £Qk/(2r),'which is determined by the spectrum of the 
forcing (Qk is the spectral power of the forcing covariance, cf. Appendix A). 
However, this is not true in NL which includes eddy-eddy interactions producing 
enstrophy/energy cascades. For example, in NL an isotropic ring forcing is spread 
as time progresses, becoming concentrated at lower wavenumbers and forming the 
characteristic dumbbell shape seen in /3-plane turbulence simulations (cf. Vallis 
and Maltrud (1993)) and consequently the homogeneous turbulent state is no 
longer characterized by the spectrum of the forcing. We can take account of this 
modification of the spectrum by performing S3T stability on the homogeneous 
state under the equivalent forcing covariance, 

« L = |(li| 2 ), (5.7) 

which maintains the observed NL spectrum, ^|Ck| 2 ^ hi the S3T dynamics. The 
NL modified eddy vorticity spectrum, ^|Ck| 2 ^b is obtained from an ensemble of 
NL simulations. Plots of ^|(k| 2 ^), under IRFn forcing are shown in Figs. 5.5a-c 
for various energy input rates, £, and damping rates, r. The departure of the 
NL spectra from the spectra of the QL and S3T equilibria is evident and this 
departure depends on the amplitude of the forcing, £, and the damping, r. 

We now demonstrate that while the fundamental qualitative prediction of 
S3T that jets form as a bifurcation in turbulence forcing and in the absence 
of turbulent cascades is verified in both QL and NL, a necessary condition for 
obtaining quantitative agreement between NL and both S3T and QL dynamics 
is that the equilibrium spectrum used in the S3T and QL dynamics be close 
to the equilibrium spectrum obtained in NL so that the stability analysis is 
performed on similar states. In the case with IRFn and r = 0.01, formation of 
persistent finite amplitude zonal jets occurs in the NL simulations at £ = 2.8£ C5Z 
(cf. Fig. 5.2b). In agreement, S3T stability analysis on the NL modified equilibrium 
IRFn spectrum (denoted S3Tb and shown in Fig. 5.5c) predicts instability for 
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Figure 5.5: Panels (a)-(d): Equilibrium enstrophy spectrum, log ((|Ck| 2 )), of NL simulations, 
in which eddy-eddy interactions are included and the k x — 0 component is excluded, for various 
damping rates, r. The example is for IRFn forcing at e = 2e CjZ . Shown are spectra for: (a) r = 1, 
(b) r = 0.1 and (c) r = 0.01. The critical e CjZ is a function of r and is obtained from S3T for each 
value of r. All spectra have been normalized. The equilibrium spectrum of the S3T (identical 
to QL) is shown in panel (d). This figure shows that for strong damping the spectrum in NL 
simulations is close to the S3T spectrum while for weak damping the equilibrium spectrum in 
NL differs substantially from that in S3T. In all cases [3 = 10. Panel (e): S3T growth rates, s r , 
as a function of the meridional wavenumber, n y , for the nonlinearly modified spectrum shown 
in panel (c) (r = 0.01). Shown are cases for £ = 2e: CjZ , e = 2.8£ c ,z and e = 10£ CjZ . It can be 
seen that S3T stability analysis forced by this spectrum predicts that jets should emerge at 
e — 2.8£ C;Z with n y = 6. S3T predictions are verified in NL as shown in the bifurcation diagram 
in Fig. 5.2b (denoted as S3Tb). 
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e > 2.8£ C5Z (cf. Fig. 5.5e). Moreover, S3T stability analysis with the S3Tb 
spectrum predicts jet formation at n y — 6 and in agreement with this prediction 
jets emerge in NL with n y = 6. Hovmoller diagrams demonstrating similar jet 
evolution in NL under IRFn forcing and in S3T under S3Tb forcing are shown in 
Fig. 5.6. We also note that agreement between NL and S3T in predictions of jet 
amplitude at large supercriticality is also obtained by using the S3Tb spectrum 
(cf. Figs. 5.2a and 5.2b). 2 



Figure 5.6: Hovmoller diagrams of U(y,t) comparing jet emergence and equilibration in an 
NL simulation under IRFn forcing (panel (a)) with an S3T simulation under S3Tb forcing (panel 
(b)). The corresponding time mean jets are shown in panel (c). This figure shows that the 
S3Tb modification of the forcing spectrum suffices to obtain agreement with NL. Parameters are 
£ — 10e c ,z, P = 10, r — 0.01. 


This influence of the eddy spectrum on jet dynamics is revealed in the case of 
IRFn at energy input rate e — 2e c ^^ shown in Fig. 5.7. Although at this energy 
input rate S3T under IRFn is structurally unstable, no jets emerge in NL. We have 
shown that agreement in bifurcation structure is obtained between NL and S3T 
when S3T analysis is performed with the S3Tb spectrum. We now examine the 
development of the NL spectrum towards S3Tb and demonstrate the close control 
exerted by this evolving spectrum on S3T stability. The evolving spectrum, shown 
in Fig. 5.8a-f, is obtained using an ensemble of NL simulations, each starting from 
a state of rest and evolving under a different forcing realization. A sequence of 

2 The spectral peaks near the k y axis do not directly influence the stability of the NL modified 
spectrum, which is determined by the distorted and broadened ring spectrum. However, while 
the spectral peaks do not influence the stability directly, they do influence it indirectly by 
distorting the incoherent spectrum. 
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S3T stability analyses performed on this evolving ensemble spectrum is show in 
Fig. 5.8g. The weak NL ensemble spectrum at t — 1 does not support instability, 
but by t = 20 the ensemble spectrum, having assumed the isotropic ring structure 
of the forcing, becomes S3T unstable. This structural instability results in the 
formation of an incipient n y = 6 jet structure which is evident by t = 50 in 
the NL simulation shown in Fig. 5.7. As the spectrum further evolves, the S3T 
growth rates decrease and no jet structure is unstable for t > 120, and decay rates 
continue to increase until t — 250 (cf. Fig. 5.8g). This example demonstrates the 
tight control on S3T stability exerted by the spectrum. Furthermore, it shows 
the close association between S3T instability and the emergence of jet structure 
in NL. 

5.5 Influence of non-zonal structures predicted by S3T on the 

TURBULENCE SPECTRUM AND ON JET DYNAMICS 

Despite S3T supercriticality, no persistent jets emerge in NL simulations with 
IRFn forcing in the interval £ CjZ < e < 2.8£ CjZ (cf. Fig. 5.2a). Comparisons of NL, 
QL and S3T simulations with IRFn forcing at e — 2e c ^ z are shown in Fig. 5.7. 
Instead of zonal jets, in the NL simulation prominent non-zonal structures are 
seen to propagate westward at the Rossby wave phase speed. These non-zonal 
structures are also evident in the concentration of power in the enstrophy spectrum 
at (\k x \,\ky\) = (1,7) (cf. top panels of Fig. 5.9). At this forcing amplitude these 
structures are essentially linear Rossby waves which, if stochastically forced, 
would be coherent only over the dissipation time scale 1/r. Coherence on the 
dissipation time scale is observed in the subdominant part of the spectrum as 
seen in the case of the (3,6) structure in Fig. 5.10c. However, the dominant (1,7) 
structure remains coherent over time periods far exceeding the dissipation time 
scale (cf. Hovmoller diagram Fig. 5.10b). This case represents a regime in which 
the flow is dominated by a single non-zonal structure. Both the concentration of 
power in and the coherence of this structure will be addressed below. 

When the forcing is increased to s = lO^z, a (0,6) jet structure emerges, 
suppresses the non-zonal (1,7) structure, and becomes the dominant structure. A 
prominent phase coherent non-zonal (1,5) structure propagating with the Rossby 
wave speed is also present, as shown in Fig. 5.11. A similar regime of coexisting 
jets and non-zonal structures is also evident at higher supercriticalities. An 
example is the case of the equilibrium state at e = 100£ CjZ (cf. Fig. 5.3) in which 
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Figure 5.7: Hovmoller diagrams of jet emergence in NL, QL and S3T simulations with IRFn 
forcing at e = 2e CjZ . Shown is U(y,t) for the NL (panel (a)), QL (panel (c)) and S3T (panel (e)) 
simulations and characteristic snapshots of streamfunction fields at t — 2000 for the NL and 
QL simulations (panels (b) and (d)). Notice that in the U(y,t) diagram for NL the color axis 
is scaled differently. Also shown are the equilibrium jets in the NL (dash-dot), QL (dashed), 
and S3T (solid) simulation (panel (f)). At e = 2s CjZ in the NL simulation no jets emerge but 
accumulation of energy in non-zonal structures with zonal wavenumber k x — 1 and meridional 
wavenumber k y — 7 is discernible. Parameters are /3 = 10, r = 0.01. 
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Figure 5.8: Panels (a)-(f): Evolution of the ensemble average enstrophy spectrum, ^|Ck| 2 ), for 
NL with IRFn forcing at e — 2s CjZ . Panel (g): Growth rates, s r , as a function of jet meridional 
wavenumber, n y , predicted by S3T stability analysis performed on the instantaneous spectrum 
at the times indicated in panels (a)-(f). The evolving spectrum renders the NL simulation S3T 
unstable at t ~ 20 and stabilizes it again at t ~ 120. Parameters are /3 = 10, r = 0.01. 


the energy of the flow is shared between the (0, 3) jet and the (1, 3) structure, as 
shown in Fig. 5.11. At this forcing level the (1,3) structure is not phase coherent, 
but its phase speed is still given by the Rossby wave speed. At even higher 
forcing similar non-zonal structures, referred to as zonons, have been reported to 
coexist with zonal jets while propagating phase incoherently at speeds that differ 
substantially from the Rossby wave speed (Sukoriansky et ah, 2008). These cases 
provide examples of the regime in which jets and non-zonal structures coexist. 

In order to study the dynamics of non-zonal structures within the framework of 
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Figure 5.9: The statistical equilibrium enstrophy spectrum, log ((|Ck| 2 )), for NL and QL 
simulations under IRFn forcing at e = 2 s CjZ (panels (a) and (b)) and e = 10e: c , z (panels (c) 
and (d)). For e — 2s CjZ the NL simulations do not support zonal jets and energy is seen to 
accumulate in the non-zonal structure (\k x \, \k y \) = (1,7). At e = 10e c ,z, persistent zonal jets 
emerge (cf. Fig. 5.6) suppressing the power in the non-zonal structures. Parameters: ft — 10, 
r — 0.01. 


S3T the interpretation of the ensemble mean in the S3T formulation is required: 
instead of interpreting the ensemble means as zonal means, interpret them rather 
as Reynolds averages over an intermediate time scale, cf. (2.13). As we have 
seen in chapter 3, the S3T stability of the homogeneous equilibrium state using 
this interpretation reveals that when the energy input rate reaches the value £ CjZ , 
which is the S3T stability threshold for the emergence of zonal jets, the state may 
already be unstable to non-zonal structures (cf. Fig. 3.2c,e). This can be also 
seen in the stability analysis shown in Fig. 5.12. In agreement with this stability 
analysis, the spectrum of the NL simulation shows concentration of power in these 
most S3T unstable wavenumbers (cf. Fig. 5.9). 

The dominance and persistence of the structures seen in these NL simulations 
can be understood from this stability analysis and its extension into the nonlinear 
regime. Because the stochastic forcing is white in time, the energy injection rate is 
fixed and state independent and, assuming linear damping at rate r dominates the 
dissipation, the total flow energy assumes the fixed and state independent mean 
value E m + E p = e/(2r). At finite amplitude the set of S3T unstable structures 
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Figure 5.10: Hovmoller diagrams of the non-zonal structures supported in the NL simulation 
of Fig. 5.7. Panel (a): evolution of the total perturbation streamfunction, ip(x,y = yo,t), 
at latitude yo — 7t/4. Panel (b): evolution of the dominant (\k x \,\k y \) — (1,7) structure of 
i!)(x,y = yo,t) at latitude yo = 7r/4. Almost half of the energy input to the system is captured 
and dissipated by this mode, which is phase coherent and propagates at the Rossby wave speed 
indicated by the dashed line. Panel (c): evolution of the (|/c x |, | k y \) — (3,6) structure at the 
same latitude. While this structure propagates at the Rossby wave speed it is not phase coherent. 
Parameters: IRFn forcing at £ = 2s CjZ , f3 = 10, r = 0.01. 
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Figure 5.11: Panel (a): Evolution of the mean flow energy, F m , which is concentrated at 
(0,6), the total eddy energy, E p , and the energy of the (1,5), (1,6) and (1,7) structures for the 
NL simulation with IRFn forcing at £ = 10£ c ,z, shown in Fig. 5.6. Panel (b): Evolution of the 
mean flow energy, F m , the total eddy energy, E p , as well as the energy of the (1,3), (1,5) and 
(1,6) structures for the NL simulation with IRFn forcing at £ = 100£ CjZ , shown in Fig. 5.3. The 
mean flow energy is concentrated at (0,3). In both panels the evolution of the energies is shown 
after statistical steady state has been reached. 
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Figure 5.12: Growth rate, s r , of the S3T non-zonal eigenfunction, e in x , as a function of 
zonal wavenumber n x and meridional wavenumber, n y for IRFn at e = 0.75£ c ,z (panel (a)) and 
e — 2s CjZ (panel (b)). The values at the axis, (0,n y ), give the growth rate of the corresponding 
jet perturbation. For £ = 0.75e c ,z the n x — 0 jet eigenfunctions are stable but the non-zonal 
perturbations are unstable with maximum instability occurring at n = (2,8). For e = 2£ CjZ the 
n x — 0 perturbations are unstable but the non-zonal perturbations are more strongly unstable, 
with maximum growth at n = (2, 8) and n = (1, 7). An NL simulation at £ = 2£ C;Z accumulates 
energy at (\k x \, \k y \) = (1,7) (cf. Fig. 5.9) while the vorticity field shows some accumulation at 
(\k x l | k y |) = (2,8) (cf. Fig. 5.8f). The stability boundary (s r = 0) is marked with thick solid line. 
For both panels f3 — 10 and r — 0.01. 
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equilibrate to allocate among themselves most of this energy which results in 
the dominance of a small subset of these structures. However, we find that in 
this competition a specific zonal jet structure has primacy so that even if this 
structure is not the most linearly unstable it emerges as the dominant structure. 

An attractive means for exploring the dynamics of the interaction between 
jets and non-zonal structures is changing the jet damping rate in the mean 
flow equation (5.1a) from r to r m and allowing it to assume values different 
from the perturbation damping rate, r, in (5.1b). (The same change of r to 
r m is also done in the QL system (5.3) and the S3T system (5.4).) In this way 
we can control the relative stability of jets and non-zonal structures as well as 
the finite equilibrium amplitude reached by the jet. This asymmetric damping 
may be regarded as a model for approximating jet dynamics in a baroclinic 
flow in which the upper level jet is lightly damped, while the active baroclinic 
turbulence generating scales are strongly Ekman damped. This asymmetry in the 
damping between upper and lower levels contributes to making jets in baroclinic 
turbulence generally stronger than jets in barotropic turbulence (Farrell and 
Ioannou, 2007, 2008). By appropriate choice of r and r m a regime can be obtained 
in which the zonal jet instability appears first as e increases. Because once jets 
are unstable they dominate non-zonal structures, in this regime zonal jets are the 
dominant coherent structure and S3T analysis based on the zonal interpretation 
of the ensemble mean produces very good agreement with NL. For example, a 
comparison of bifurcation structures among S3T, QL and NL under NIF and 
IRFn forcing using the asymmetric damping r — 0.1 and r m = 0.01 demonstrates 
that jets emerge at the same critical value in S3T, QL and NL (cf. Figs. 5.13a and 
5.13b). This agreement, which has been obtained by suppression of the non-zonal 
instability up to £ C5Z , implies that in the simulations with symmetric damping the 
disagreement in the S3T prediction for the first emergence of jets (cf. Fig. 5.2) can 
be attributed to modification of the background spectrum by the prior emergence 
of the non-zonal structures. Moreover, zonal structures once unstable immediately 
dominate non-zonal structures assuring that S3T dynamics based on the zonal 
mean interpretation of the ensemble mean produces accurate results. 

A comparison of the development of jets in S3T, QL, and NL with this asym¬ 
metric damping and NIF forcing, shown in Figs. 5.14 and 5.15, demonstrates 
the accuracy of the S3T predictions. S3T stability analysis predicts that in 
this case with NIF forcing maximum instability occurs at n y — 6. When these 
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maximally growing eigenfunctions are introduced in the S3T system the jets grow 
exponentially at first at the predicted rate and then equilibrate. Corresponding 
simulations with the QL and NL dynamics reveal nearly identical jet growth 
followed by finite amplitude equilibration (shown in both Figs. 5.14 and 5.15). 
Similar results are obtained with IRFn forcing. This demonstrates that the S3T 
dynamics comprises both the jet instability mechanism and the mechanism of 
finite amplitude equilibration. 

Although no theoretical prediction of this bifurcation behavior can be made 
directly from NL or QL, they both reveal the bifurcation structure obtained from 
the S3T analysis. By suppressing the peripheral complexity of non-zonal structure 
formation by non-zonal S3T instabilities, these simulations allow construction 
of a simple model example that provides compelling evidence for identifying jet 
formation and equilibration in NL with the S3T theoretical framework. Moreover, 
agreement among the NL, QL and S3T bifurcation diagrams shown in Figs. 5.13a 
and 5.13b provides convincing evidence that turbulent cascades, which are not 
present in S3T or QL, are not required for jet formation. 

While under NIF agreement between NL and S3T equilibrium jet amplitudes 
extends to all values of £, under IRFn the NL and S3T equilibrium amplitudes 
diverge at larger values of e (cf. Figs. 5.13a and 5.13b). This difference among 
NL, QL and S3T at large e cannot be attributed to nonlinear modification of the 
spectrum, which is accounted for by use of the S3Tb spectrum (cf. S3Tb response 
in Fig. 5.13b). Rather, this difference is primarily due to nonlinear eddy-eddy 
interactions retained in NL that disrupt the up-gradient momentum transfer. This 
disruption is accentuated by the peculiar efficiency with which the narrow ring 
forcing, IRFn, gives rise to vortices, as can seen in Fig. 5.1d-f. The more physical 
distributed forcing structures do not share this property (cf. Fig. 5.1). We verify 
that the narrow ring IRFn forcing is responsible for depressing NL equilibrium jet 
strength at high supercriticality by broadening the forcing distribution to assume 
the form IRFw (cf. Appendix H as well as Fig. 5.1 for IRFn-IRFw comparison). 
Using IRFw while retaining other parameters as in Fig. 5.13b, we obtain agreement 
between S3T, QL and NL simulations, as is shown in Fig. 5.13c. 
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5.6 Identification of intermittent jets with stable S3T zonal eigen¬ 
functions 

For subcritical forcing S3T predicts a stable homogeneous statistical equilibrium 
and a set of eigenfunctions that govern the decay of perturbations to this equilib¬ 
rium. We wish to show that these eigenfunctions are excited in NL by fluctuations 
in the turbulence and that this excitation gives rise in NL simulations to the 
formation of intermittent jets with the form of these eigenfunctions. 

As an example, consider the simulation with asymmetric damping and IRFn 
subcritical forcing shown in Fig. 5.16. For these parameters the least damped 
eigenfunctions are zonal jets and confirmation that the intermittent jets in NL, 
shown in the top panel of Fig. 5.16, are consistent with turbulence fluctuations 
exciting the S3T damped modes is given in the bottom panel of Fig. 5.16 where the 
intermittent jets resulting from stochastic forcing of the S3T modes themselves are 
shown. This diagram was obtained by plotting U(y,t) — Re a n y (t)e iriyy , 

with OL Uy independent red noise processes, associated with the damping rates, 
|s(n)|, of the first N = 15 least damped S3T modes. These a Uy are obtained 
from the Langevin equation, da ny /dt = s(n y ) a Uy + £(£), with £(t) a 5-correlated 
complex valued random variable. 

The fluctuation-free S3T simulations reveal persistent jet structure only coinci¬ 
dent with the inception of the S3T instability, which occurs only for supercritical 
forcing. However, in QL and NL simulations fluctuations excite the damped 
manifold of modes predicted by the S3T analysis to exist at subcritical forcing 
amplitudes. This observation confirms the reality of the manifold of S3T stable 
modes. 

In NL and QL simulations these stable modes predicted by S3T are increasingly 
excited as the critical bifurcation point in parameter space is approached, because 
their damping rate vanishes at the bifurcation. The associated increase in zonal 
mean flow energy on approach to the bifurcation point obscures the exact location 
of the bifurcation point in NL and QL simulations compared to the fluctuation-free 
S3T simulations for which the bifurcation is exactly coincident with the inception 
of the S3T instability (i.e. Fig. 5.13a, 5.13b and 5.13c). 
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5.7 Verification in NL of the multiple jet equilibria predicted by 
S3T 

As is commonly found in nonlinear systems, the finite amplitude equilibria 
predicted by S3T are not necessarily unique and multiple equilibria can occur 
for the same parameters. S3T provides a theoretical framework for studying 
these multiple equilibria, their stability and bifurcation structure. An example 
of two such S3T equilibria are shown in Fig. 5.17 together with their associated 
NL simulations. As the parameters change these equilibria may cease to exist 
or become S3T unstable. Similar multiple equilibria have been found in S3T 
studies of barotropic /3-plane turbulence (Farrell and Ioannou, 2003, 2007; Parker 
and Krommes, 2014a) and in S3T studies of baroclinic turbulence (Farrell and 
Ioannou, 2008, 2009c) and the hypothesis has been advanced that the existence 
of such multiple jet equilibria may underlie the abrupt transitions found in the 
record of Earth’s climate (Farrell and Ioannou, 2003; Wunsch, 2003). 

5.8 Conclusions 

In this chapter predictions of S3T for jet formation and equilibration in barotropic 
/3-plane turbulence were critically compared with results obtained using QL and 
NL simulations. Throughout this chapter the zonal mean-eddy decomposition 
(section 2.2) was used for all three NL, QL and S3T systems. The qualitative bifur¬ 
cation structure predicted by S3T for emergence of zonal jets from a homogeneous 
turbulent state was confirmed by both the QL and NL simulations. Moreover, 
the finite amplitude equilibrium jets in NL and QL simulations were found to 
be as predicted by the fixed point solutions of S3T. Differences in jet formation 
bifurcation parameter values between NL and QL/S3T were reconciled by taking 
account of the fact that the spectrum of turbulence is substantially modified in 
NL. Remarkably, the modification of the spectrum in NL could be traced in large 
part to emergence of non-zonal structures through S3T instability. When account 
is taken of the modification of the turbulent spectrum resulting substantially from 
these non-zonal structures, S3T also provides quantitative agreement with the 
threshold values for the emergence of jets in NL. The influence of the background 
eddy spectrum on the S3T dynamics was found to be immediate, in the sense that 
in spin-up simulations jets emerge in accordance with the instability calculated 
on the temporally developing spectrum. The fact that jets are prominent in 
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observations is consistent with the robust result that when a jet structure emerges 
it has primacy over the non-zonal structures, so that even if the jet eigenfunction 
is not the most linearly S3T unstable eigenfunction, the jet still emerges at finite 
amplitude as the dominant structure. 

These results confirm that jet emergence and equilibration in barotropic /3-plane 
turbulence results from the cooperative quasi-linear mean flow-eddy instability 
that is predicted by S3T. These results also establish that turbulent cascades 
are not required for the formation of zonal jets in /3-plane turbulence. Moreover, 
the physical reality of the manifold of stable modes arising from cooperative 
interaction between incoherent turbulence and coherent jets, which is predicted 
by S3T, was verified in this work by relating observations of intermittent jets in 
NL and QL to stochastic excitation by the turbulence of this manifold of stable 
S3T modes. 

5.9 Bibliographical note 

This chapter is an adaptation from the paper by Constantinou et al. (2014a). The 
NIF forcing used in this chapter was first used by Williams (1978) in order to to 
parametrize excitation of barotropic dynamics by baroclinic instabilities. It was 
also used by DelSole (2001) in his study of upper-level tropospheric jet dynamics 
and in the study of jet formation using S3T dynamics by Farrell and Ioannou 
(2003, 2007) and Bakas and Ioannou (2011). The isotropic narrow ring forcing, 
IRFn, has been used extensively in studies of /3-plane turbulence (cf. Vallis and 
Maltrud, 1993) and was also used in the recent study of Srinivasan and Young 
(2012). It was introduced by Lilly (1969), in order to isolate the inverse cascade 
from the forcing in a study of two dimensional turbulence. 
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Figure 5.13: Bifurcation structure comparison for jet formation in S3T, QL, and NL with 
asymmetric damping. Shown is the zmf index of jet equilibria for NIF (panel (a)), IRFn 
(panel (b)) and IRFw (panel (c)) as a function of the forcing amplitude s/£ c ,z for the NL 
simulation (dash-dot and circles), the QL simulation (dashed and dots) and the corresponding 
S3Ta simulation (solid and diamonds). Also shown in panel (b) is the zmf that is obtained from 
S3T simulations forced with the nonlinearly modified S3Tb spectrum (calculated from ensemble 
NL simulations at £ = 20e c ,z)- Parameters are (3 = 10, r — 0.1, r m = 0.01. 
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Figure 5.14: Hovmoller diagrams of jet emergence in the NL, QL and S3T simulations for 
NIF at £ = 1.5e c ,z with asymmetric damping. Shown is U(y,t) for the NL (panel (a)), QL (panel 
(c)) and S3T (panel (e)) simulations and also characteristic snapshot of the vorticity fields at 
t = 2000 for NL and QL simulations (panels (b) and (d)). Also shown are the equilibrium jets in 
the NL (dash-dot), QL (dashed), and S3T (solid) simulation (panel (f)). This figure shows that 
S3T predicts the structure, growth and equilibration of weakly forced jets in both the QL and 
NL simulations. Parameters are: [3 — 10, r — 0.1, r m = 0.01. 
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Figure 5.15: Same as Fig. 5.14 but with forcing at £ = 20£ c ,z- While initially jets emerge 
having the structure of the most unstable S3T mean flow eigenfunction with n y =6, at later 
times, following a series of mergers they equilibrate to a finite amplitude state with n y = 4. 
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Figure 5.16: Hovmoller diagrams of intermittent jet structure in NL and QL simulations 
at subcritical forcing s — 0.8s c ,z- Shown are U(y,t) for NL (panel (a)) and QL (panel (b)) 
simulations and the U(y,t) that results from random excitation of the S3T damped modes (panel 
(c)). These plots were obtained using IRFn forcing with r — 0.1, r m = 0.01. This figure shows 
that the manifold of S3T damped modes are revealed by being excited in the fluctuating NL 
and QL simulations. Planetary vorticity gradient: /5 = 10. 
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Figure 5.17: Realizations in NL simulations of multiple equilibria predicted by S3T. Show 
are Hovmoller diagrams of NL simulations showing the equilibrium with 4 jets (panel (a)) and 
with 5 jets (panel (c)). Also shown is comparison of the S3T equilibrium jets (solid) with the 
average jets obtained from the NL simulation (dashed) for the two equilibria (panels (b) and 
(d)). Parameters: NIF forcing at e = 10£ c ,z, r = 0.1, r m = 0.01 and /3 — 10. 
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6 

S3T stability of inhomogeneous turbulent 

equilibria 


We have already seen that with homogeneous forcing the S3T system has homoge¬ 
neous equilibria for any level of forcing, but also inhomogeneous equilibria in the 
form of zonal jets, as those obtained in chapter 5, or non-zonal inhomogeneous 
equilibria in a moving frame of reference, as the traveling wave solutions that 
were obtained in chapter 4 (cf. Fig. 4.7). In chapter 3 we presented systematic 
methods that enabled us to determine the stability of the homogeneous state and 
we were able to predict the critical parameters for which the symmetry of the 
homogeneous state is broken. The stability of finite amplitude zonal jet equilib¬ 
ria to zonal jet perturbations has been already studied by Farrell and Ioannou 
(2003) and more recently by Parker and Krommes (2014a). Here we present more 
general methods for determining the stability of inhomogeneous states to zonal 
but also non-zonal perturbations. With this more general stability analysis we 
demonstrate that the phenomenon of jet merging is properly understood as an 
S3T instability and consequently this phenomenon is properly understood in the 
framework of statistical state dynamics. We also show that the transition from 
zonal to non-zonal turbulent states is also predicted by S3T stability analysis. 
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6.1 Stability of finite amplitude zonal jet S3T equilibria to zonal 

JET MEAN FLOW PERTURBATIONS 

Consider first the equilibrium states that arise in the simpler S3Tz system (2.21), 
i.e., the S3T system in which ensemble means are interpreted as zonal means. In 
S3Tz by construction the mean flows are zonal jets and the associated equilibria, 
when they exist, can only be zonal jets. The zonal jet equilibria arise as a 
bifurcation of the homogeneous equilibrium state that becomes S3T unstable for 
energy injection rates, £, that exceed a critical value, £ c?z (see Figs. 5.2 and 5.13). 
The resulting zonal jet equilibria have the characteristic property that the number 
of jets decreases as the supercriticality, s/£ C: z, increases 1 * ; examples are shown in 
Figs. 5.3, 5.6, 5.14 and 5.15. We wish to study the stability of these jet equilibria in 
order to understand the mechanism underlying the transition from one equilibrium 
state to another. In order to proceed with the stability analysis of the equilibria, 
we must first determine the equilibrium solutions, (U e (y), C e (x a — x b , y a , y^)), 
with adequate accuracy in order to obtain good estimates of their stability. While 
stable equilibrium solutions can be in principle obtained with the required accuracy 
with forward time-integration of the S3T system, forward time-integration cannot 
determine unstable equilibria and consequently we must resort to continuation 
methods in order to obtain all the fixed points of the S3T equations. Both stable 
and unstable equilibrium states can be determined with great accuracy and ease 
using the continuation methods described in Appendix I. 

As discussed in chapter 2, the linear stability of S3T equilibria is studied 
through eigenanalysis of the operator governing the linearized S3T evolution of 
the perturbations (5Z, 8C ) about the equilibrium state: 

dt 8Z = A e 8 Z + K{8C) , (6. la) 

dt 8C ab — (yia + a£) 8C ab + ( '8A a + SAbj C^ b , (6.1b) 

with A e = A(U e ) and 8A = A(\J e + 8JJ) — A e . Note that discretized with N 
points the dimension of the perturbation state in (6.1) is 0(7V 4 ). For example, 
if we use a modest discretization grid of N x = N y = 2 6 , the dimension of the 
equivalent matrix operator governing the stability of (6.1) is 2 24 x 2 24 ~ 10 7 x 10 7 . 
Despite the enormity of the size of the operators the real part, s r , of the maximally 

1 £ c ,z is the minimum energy input rate for the instability of the homogeneous state to zonal 

jets. 
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growing eigenvalue, s, and the corresponding spatial structure of the eigenfunction, 
(SZ, SC), can be still obtained numerically using the power method. The imaginary 
part of the eigenvalue, si, can then be determined by solving (6.1a): 

(s r + isi + A e ) SZ = K(6C) , (6.2) 

for Si. This procedure is still computationally expensive because it requires 
time-integration of a state vector of dimension N x N y + N^Ny. The dimension 
of the system can be reduced by a square root when the equilibrium states are 
zonal, as Bloch’s theorem (cf. Appendix D) requires that the spatial structure of 
the eigenfunction can be assumed to be of the form: 

8Z(x) = e mxX 8Z nx (y) , (6.3a) 

8C(x a , x b ) = e inx ( Xa+Xb ' ) / 2 8C nx (x a - x b , y a , y b ) , (6.3b) 

and the stability of the equilibrium is determined by evolving each zonal wavenum¬ 
ber n x separately. More details regarding the method for determining the stability 
of (6.1) are discussed in Appendix J. While for the case of the homogeneous 
equilibria the eigenfunctions are single harmonics, i.e., 5Z Ux (y ) = e inyV (cf. (E.4)), 
for zonal jet S3T equilibria the eigenfunctions are single harmonic only in x and 
have full spectrum in y (cf. Fig. 6.7 and also Fig. 6.12). However, Bloch’s theorem 
(cf. Appendix D) restricts the meridional structure of the eigenfunction 6.3. For 
example, the mean flow perturbation must be of the form 

SZ nx (y) = e iqvV g(y) , (6.4) 

with g(y) any function that has the same periodicity as the equilibrium jet, U e (y ), 
and similarly for the perturbation covariance eigenfunction, SC Ux (cf. (D.7)). For 
a n y -jet equilibrium function g is of the form: g(y) = e irnnyV . Wavenumber q y 
is called “Bloch wavenumber” and takes all values q y < n y / 2 (for our channel of 
length L x — L y = 2n wavenumber q y takes all integer values q y < n y /2). Therefore 
a, q y — 0 Bloch eigenfunction will have power at wavenumbers 0, ±n y , ±2 n y , ..., 
while ag y = l Bloch eigenfunction will have power at wavenumbers ±1, ±{n y ± 
1),±(2 % ±1),.... 

Using the continuation methods described in Appendix I we find a series of 
zonal jet equilibria that are characterized by different number of prograde jets, 
n y . Consider for example the case with n y — 6 jets and for the parameters: 
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NIF forcing, /? = 10, r = 0.1 and r m = 0.01 (cf. chapter 5, case presented in 
Fig. 5.13a). For these parameter values the homogeneous equilibrium becomes 
first unstable to zonal jets with n y — 6 at s = £ c , z and inhomogeneous zonal 
jet equilibria exist for all energy input rates e > £ CjZ that do not exceed 686£ c?z . 
The nonexistence of zonal equilibria with n y = 6 for s > 686£ C5Z is attributed 
to the inability to find equilibria when the flow starts supporting stable modal 
structures that produce strong vorticity fluxes in the neighborhood of their critical 
layers. Recall (cf. Appendix B) that although all the S3T equilibria are necessarily 
hydrodynamically stable they may be S3T unstable. Specifically, the n y — 6 
equilibria become S3T unstable for e > 20s CjZ , that is for energy input rates 
substantially lower than the energy input rate at which the equilibria cease to 
exist. Some stable and unstable S3T equilibria, their associated planetary vorticity 
gradient, (3 — Uy y , as well as the amplitude of the jets as a function of energy input 
rate are shown in Fig. 6.1. Note that because of the presence of dissipation, the 
flow can remain stable although the mean planetary vorticity gradient changes 
sign and becomes slightly negative in limited regions when the jet is retrograde 
(cf. Fig. 6.1c). Note also that as the energy input rate is increased the amplitude 
of the mean flow grows very gradually after e/e c , z ~ 100 (cf. Fig. 6.Id) and the 
energy input of the stochastic forcing is absorbed mainly by the perturbation field, 
which indicates that the extra energy in the perturbation field is not communicated 
to the mean flow. This happens because for large enough energy input rates (for 
e > 20e CjZ for n y — 6) a nearly neutral modal structure is supported by the flow, 
with critical layers where (3 — U yy vanishes, and this mode absorbs most of the 
incoming energy without producing appreciable upgradient vorticity fluxes to 
support a stronger flow. This modal structure eventually develops strong critical 
layers and strongly localized vorticity fluxes that make the existence of an S3T 
equilibrium impossible. 

A comprehensive mapping of S3T zonal jet equilibria together with their S3T 
stability as a function of the supercriticality, £/£ C;Z , and the number of jets, is 
shown in the balloon diagram Fig. 6.2. Jet equilibria exist in the yellow region of 
the diagram. For values of £/£ c , z and n y below the lower bounding curve (dashed) 
the only S3T equilibrium is the homogeneous state with no mean flow and the 
dashed line is the curve of neutral S3T stability of the homogeneous state. For 
values of e and n y above the upper bounding curve no jet equilibria exist. The 
S3T stability of certain jet equilibria to zonal jet perturbations (i.e. with n x — 0 
in (6.3)) is indicated with a closed circle when it is stable and with an open circle 
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Figure 6.1: (a) S3T zonal flow equilibria with n y = 6 jets for supercriticality s/e c ,z — 

3,10, 20, 50, 200. The equilibria for s/e c , z = 20, 50 (dashed) are S3T unstable, and those with 
e/e c ,z < 20 are S3T stable, (b) The potential vorticity gradient, [3 — U yy for the corresponding 
equilibria shown in panel (a), (c) The maximum growth rate of A(U e ) for the zonal flow equilibria 
with n y = 6 jets as a function of the supercriticality. All equilibria are hydrodynamically stable. 
S3T equilibria with n y = 6 jets cease to exist at e/e c ,z = 686. (d) the amplitude of the 

equilibrium jets as a function of e/e c ,z . The amplitude of the jets does not increase substantially 
for e/e CjZ > 100 and the extra energy that is imparted in the flow is absorbed by the perturbation 
field without being communicated to the mean flow. Parameters: NIF forcing, j3 — 10, r = 0.1 
and r m — 0.01. 
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Figure 6.2: S3T zonal flow equilibria as a function of the number of prograde jets, n y , 
and the forcing amplitude e/e c ,z- Jet equilibria exist in the yellow shaded region. The lower 
bounding curve is the curve of neutral stability of the homogeneous state (the homogeneous 
state is S3T stable for all e for n y > 12). For e and n y above the upper bounding dashed-dot 
curve jet equilibria do not exist. S3T stable finite amplitude jet equilibria to jet perturbations 
are indicated with a full circle, S3T unstable equilibria to jet perturbations are indicated with an 
open circle. For a range of £ there exist multiple S3T stable equilibria characterized by a different 
number of jets. Near the curve of marginal stability only the jet equilibrium that corresponds to 
the maximal instability of the homogeneous state, n y = 6 is stable, while the neighboring jet 
equilibria are Eckhaus unstable. Parameters: NIF forcing, — 10, r = 0.1 and r m = 0.01. 


The balloon diagram shows that for a range of values of e multiple stable 
equilibria exist (these correspond to multiple climate states in this barotropic 
model). As e is increased all n y -jet equilibria become eventually S3T unstable. 
When they become unstable the turbulent flow reorganizes, the mean flow merges 
and transitions to an available S3T jet stable equilibrium with fewer jets. For 
example, at e/e c ^ = 10 the equilibria with n y = 3,4, 5, 6 jets are all stable 
(in chapter 5 we have seen that this is also reflected in the NL simulations; 
cf. Fig. 5.17). If we increase the supercriticality to a value e /£ c ,z > 10 the n y = 6 
jet becomes S3T unstable and the turbulent flow reorganizes to a state with 
n y — 3,4, 5 or even n y — 2 jets. This implies that if the energy input rate were 
to increase in a 6-jet equilibrium the turbulent state would transition to a state 
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with a fewer number of jets through a process of jet mergers. 

We now give an example which demonstrates in an actual S3T simulation 
that jet mergers do not occur because of the hydrodynamic instability of the 
jets but rather due to their S3T instability. In the simulation shown in Fig. 6.3, 
for e/e c , z — 100. The flow evolves towards jet configurations with fewer and 
stronger jets though a sequence of jet mergers. Interestingly, the jet states of 
the flow when they do not merge evolve slowly staying close to corresponding 
S3T unstable equilibria till finally they get attracted to the first available S3T 
stable equilibrium (for the chosen parameters the first S3T stable equilibrium 
to jet perturbations {n x = 0) has n y — 3, all equilibria with n y > 3 are S3T 
unstable to n x = 0 perturbations; cf. Fig. 6.2). Because the evolution of the jets 
when they do not merge is slow we can interpret the jet merging process as an 
instability of the time evolving state of the system and because the strength of 
the hydrodynamic instability of a barotropic flow depends on the strength of 
the violation of the Rayleigh-Kuo criterion (3 — U yy , it is natural to attribute jet 
merging to the hydrodynamic instability of the flow. We see also, that during 
this evolution the flow stays close to the various S3T equilibria before moving 
to the next equilibrium. The maximum growth rate of A(U) that governs the 
hydrodynamic stability of the flow is negative at each time indicating that jet 
merging can not be attributed to the hydrodynamic stability of the flow. Thus 
we conclude that the flow reorganization from a state with 4 jets to the state 
with 3 jets does not occur due to hydrodynamic instability of the zonal flow but 
is an inherent S3T phenomenon. In the specific example because the flow stays 
close to the various S3T unstable equilibria jet merging can be attributed to the 
S3T instability of these unstable equilibria that behave as unstable saddles of the 
evolving flow. 

As we increase the energy input rate, £, the stable S3T equilibria have fewer 
jets and the structure of the equilibrium zonal flow equilibria, t/ e , acquires a 
particular shape; see Fig. 6.4. While just above the stability boundary the jets 
are to a good approximation sinusoidal, U e ~ sin (n y y) (cf. Fig. 6.4a), at higher 
£ the retrograde parts of the jets become parabolic while the prograde parts of 
the jets become increasingly pointed. This particular structure closely resembles 
the shape of the observed jets in planetary atmospheres (see the 24°N jet on 
Jupiter, shown in Fig. 1.4c as well as the parabolic equatorial jets in Uranus and 
Neptune, shown in Figs. 1.14a,b). That the retrograde parts of the jets are nearly 
parabolic is in agreement with potential vorticity (PV) mixing or homogenization 
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£ = 100s 



Figure 6.3: (a) Hovmoller diagram of the zonal mean flow U(y,t) showing a series of jet 

mergers leading finally to a state with 3 jets. Solid line marks the zero contour. The S3T 
simulation starts from the homogeneous equilibrium state perturbed by a random mean flow 
perturbation, (b) The corresponding evolution of the mean flow energy, E m , and perturbation 
energy, E p . Marked also are the energies of the homogeneous equilibrium and of the jet equilibria 
with n y = 2,3,4. The flow for long periods is close to S3T unstable equilibria, till it finally 
settles to the first available stable S3T equilibrium. For these parameters (cf. Fig. 6.2) the 
n y — 1,2 jet equilibria are also stable, but both have larger/smaller mean flow/perturbation 
energy, (c): The evolution of the maximum growth rate of A(JJ) for the instantaneous U(y, t). 
The flow is at every time instant hydrodynamically stable. Parameters are as in Fig. 6.2 with 
forcing at £ = 100£ c ,z- 


110 





































arguments (Baldwin et al., 2007; Dritschel and McIntyre, 2008; Dritschel and 
Scott, 2011; Scott and Dritschel, 2012). According to these mixing arguments 
the primary process in barotropic turbulence is the irreversible mixing by the 
strongly nonlinear processes that leads to homogenized regions of the mean PV, 
q = /3y — U y , that manifest as a staircase in a diagram of the PV as function 
of latitude (see also the discussion in section 1.2.1 and Fig. 1.13). In Fig. 6.6 
we plot the PV structure for the S3T equilibrium jets shown in Fig. 6.4. It is 
evident that S3T dynamics do produce a staircase structure when the forcing is 
strong. This shows that PV staircases are produced by S3T dynamics despite 
the absence of all eddy-eddy interactions. Note that the staircase structure 
obtained here is similar to the staircases obtained in the fully nonlinear and 
nearly inviscid simulations of Scott and Dritschel (2012) as well as those that are 
observed in experiments (cf. Fig. 6.6). These S3T equilibria with PV staircase 
structure provide a counterexample to the necessity of wave breaking and strong 
nonlinearity for the formation of the staircase structure in barotropic turbulence 
argued by McIntyre and collaborators. 

The primary mechanism responsible for the specific shape of the S3T equilibrium 
flows at strong forcing is that turbulence acts on the mean flow as negative diffusion 
(cf. section 3.4) and that the flow must be hydrodynamically stable, which means 
that the necessary Rayleigh-Kuo criterion for stability cannot be violated (in the 
limit of infinitesimal friction, i.e., £* = ek'j/r^ oo). The latter requirement 
constrains only the retrograde regions of the flow (where U yy > 0) and brings 
/? — U yy to the minimum value that does not violate Rayleigh-Kuo, i.e, brings 
/3 — U yy to zero, making in this way the flow parabolic, while the former requirement 
leads to the formation of nearly linear prograde flows joined with sharp wedge-like 
peaks. As the prograde jets become sharper, the zonal mean flow energy spectrum 
develops a slope (cf. Fig. 6.4), as is expected from the near discontinuity of 
the derivative of the prograde jets (a discontinuity would predict a k~ A zonal 
energy spectrum; see also discussion in section 1.2.1). The absence of eddy-eddy 
interactions in the S3T dynamics leads us to conclude that the observed zonal 
energy spectrum cannot be attributed to the anisotropic and incoherent inverse 
turbulent energy cascade, as it was recently proposed by Galperin et al. (2004). 

For very high supercriticalities only finite amplitude states with jets having the 
largest allowed scale, n y — 1, exist (this occurs in the balloon diagram Fig. 6.2 for 
e/e c ^ > 800). Higher supercriticalities cannot sustain S3T fixed points because 
the periodic box does not allow jets larger than the box size and the S3T dynamics 
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Figure 6.4: Panels (a,d,g,j): The structure of the S3T stable equilibrium zonal mean flow, 
f/ e , for excitation amplitudes e/sc^ — 2,8,150,800,10 4 . Panels (b,e,h,k): The corresponding 
mean vorticity gradient, ft — U yy . Dash-dotted line marks the planetary vorticity gradient in the 
absence of mean flow, ft — 10. Panels (c,f,i,l): The energy spectrum of the equilibrium zonal 
mean flow together with the k y 5 slope. For the highly supercritical jets the energy spectrum 
has approximately an k~ b dependence. It is argued that in the inviscid limit this slope should 
approach k ~ 4 as the prograde jet becomes increasingly sharp. Other parameters as in Fig. 6.2. 



Figure 6.5: The zonal mean PV, q e (y) = fty — U y , that corresponds to the equilibrium zonal 
flows shown in Fig. 6.4. Dashed lines correspond to the planetary PV, f3y. 
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Figure 6.6: (a) The zonal mean flow and (b) the corresponding PV from the Grenoble 

rotating tank experiment. The PV staircases structure closely resembles the one in Fig. 6.5. 
(Taken from Read et al. (2007). The axes were flipped so that the orientation matches Fig. 6.5.) 
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Figure 6.7: The jet equilibria with n y — 5 and n y — 7 near the marginal curve are both 
found to be Eckhaus unstable to an n y — 6 zonal jet perturbation. (a,d) Contour plots of the 
equilibrium streamfunction, T e (y), together with the zonally averaged zonal velocity U e (y) (thick 
black line) for jet equilibria with n y — 5 jets at s/s c ,z = 1-05 and n y — 7 jets at £/£ c ,z = 1.1. (b,e) 
The structure of the mean flow streamfunction of the maximally growing S3T eigenfunctions for 
the jet equilibria together with its eigenvalue a. (c,f) The energy spectrum of the corresponding 
S3T eigenfunctions as a function of the meridional wavenumber n y . While both eigenfunction 
show maximum power at n y m 6 they have also non-zero power spectrum at other wavenumbers. 
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eventually produce chaotic and non-stationary trajectories of the statistics of the 
turbulent flow. 

We return to the equilibrium states that emerge after the homogeneous equilib¬ 
rium is broken. We note that for 1 < s/e c , z <1.2, while inhomogeneous equilibria 
with n y = 5,6,7 are found, only the equilibrium with n y = 6 is found to be 
stable. The equilibria with n y — 5 and n y = 7 are both unstable with their most 
unstable eigenfunction being a zonal jet with maximum power spectral power at 
n y — 6 (cf. Fig. 6.7). This instability is the universal Eckhaus instability that 
occurs near the neutral stability boundary and attracts all finite amplitude states 
to the structure of the most unstable eigenfunction of the homogeneous state; 
in this case n y — 6. The nonlinear dynamics near the marginal curve obey a 
Ginzburg-Landau equation (cf. Parker and Krommes (2014b)). It should be noted 
that the accuracy of the approximation of the dynamics by the Ginzburg-Landau 
equation is unfortunately limited only to parameter values that are very close 
to the stability boundary, i.e., for s < 1.0l£ c ,z- Also the jet mergers that occur 
under Ginzburg-Landau dynamics are associated with the equilibration of the 
Eckhaus instabilities and are very different from the jet mergers that are seen 
in S3T simulations at higher supercriticality (cf. Parker and Krommes (2014b) 
and Fig. 6.3). 

6.2 Stability of finite amplitude zonal jet S3T equilibria to non- 

ZONAL MEAN FLOW PERTURBATIONS 

The method described in the previous section for the stability of zonal jet S3T 
equilibria to zonal jet perturbations (i.e. mean flow perturbations with n x = 0) 
will be employed now to determine the stability of the zonal jet equilibria that 
were discussed in the previous section to non-zonal mean flow perturbations, 
(i.e. mean flow perturbations with n x ^ 0). The question we want to address 
is: can a jet be stable to jet like perturbations {n x — 0) but still be unstable 
to non-zonal large scale structures? We will demonstrate that all the stable jet 
equilibria of the balloon diagram Fig. 6.2 are unstable to non-zonal perturbations. 
We note that these equilibria were obtained with NIF forcing and with the eddy 
flow being damped more strongly than the large-scale flow (both the jet and 
non-zonal large-scale components). We have obtained similar results when both 
flows were dissipated equally and the forcing was isotropic. Therefore we believe 
the results presented here are not pathological. This instability of strong jet 
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equilibria to non-zonal perturbations is intriguing because our theory at the 
level of our numerical implementation predicts that in barotropic turbulence the 
statistical equilibria must have a non-zonal component, contrary for example to 
the observations on Jupiter which strongly suggest that the equilibria are almost 
purely zonal. 

Consider first the stability of two of the n y — 6 jet equilibria shown in Fig. 6.1, 
the equilibrium for e/e c , z = 10, which is stable to jet perturbations, and the 
equilibrium for e/e^ z = 200, which is unstable to jet perturbations. The cor¬ 
responding S3T growth rates for zonal and non-zonal perturbations are shown 
in Fig. 6.8. In both cases the equilibria are unstable to non-zonal mean flow 
perturbations with 1 < n x < 7. We have also plotted the growth rates that would 
obtain if the equilibrium were homogeneous and the eddy field had the eddy 
energy zonal spectrum of the S3T equilibrium with jets (shown in Fig. 6.8c,d). 
We see that in general the presence of the jet increases the S3T stability of the 
turbulent state. However, this increase of S3T stability with jet strength does 
not eliminate the non-zonal instability of strong jets. For example consider the 
instability of the strong jet equilibria with n y — 1,2 that are stable to n x = 0 jet 
perturbations. The n y — 1 equilibrium state at e/e CjZ = 1000 is shown in Fig. 6.9 
together with the streamfunction of the least stable mode for perturbations with 
n x — 0 and the most unstable modes for perturbations with n x = 1,2. Maximum 
instability for this jet equilibrium is found at n x — 1, as shown in Fig. 6.10a and 
similarly for the n y = 2 jet equilibrium at s/e c , z = 800 in Fig. 6.10b. 

Consider now the stability of the zonal jet equilibrium that we obtained with 
non-broadband forcing in our discussion of modulational instability in chapter 4. 
In that example, the fastest growing instability of the homogeneous state is a zonal 
3-jet perturbation. In an S3T integration this zonal jet structure emerges and 
saturates to a quasi-stationary finite amplitude 3-jet flow. This flow eventually 
breaks giving way to a finite amplitude traveling wave with maximal power 
at wavenumber (n x ,n y ) = (1,4) (see Fig. 4.7). We show here that this finite 
amplitude wave with maximum power at (1,4) emerges because the 3-jet state is 
unstable to non-zonal perturbations with the property that all the instabilities 
equilibrate to a finite amplitude state with maximal power at (1,4). With the 
methods described in the previous section we first obtain the 3-jet equilibrium 
and then calculate its stability. The S3T equilibrium and the maximally growing 
eigenfunction for n x — 0,..., 6 are shown in Fig. 6.11 and the spectrum of the 
n x — 0,..., 6 eigenfunctions is shown in Fig. 6.12. The growth rate of the 3-jet 
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Figure 6.8: The growth rate of the maximally growing S3T eigenfunction of the 6-jet 
equilibrium at supercriticality e/e c ,z = 10 (panel (a)), and e/e c ,z — 200 (panel (b)) as a function 
of the zonal wavenumber n x of the perturbation (blue line, 0 )• In (a) the jet is S3T stable and 
in (b) unstable. The jet that is stable to n x = 0 becomes unstable for non-zonal perturbations. 
For comparison we also plot the growth rate of the unstable homogeneous equilibrium with 
perturbation energy equal to that of the inhomogeneous equilibrium (red line, □ ). The 
equilibrium perturbation energy associated with the 6-jet equilibrium flows is shown respectively 
in panels (c) and (d). The presence of jets is seen to generally increase the S3T stability of the 
flow. Other parameters as in Fig. 6.2. 
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Figure 6.9: Stability of the zonal jet S3T equilibrium with n y — 1 jet at e/sc^ — 1000 
(a): Contour plot of the equilibrium mean flow streamfunction, T e (y), and a plot of the zonal 
velocity, U e (y ), (thick black line). The jet maximum is 4.1 and the minimum is —3.3. The 
mean flow contains 92% of the total energy of the flow. (b,c,d): The streamfunction of the most 
unstable mode, (5T, for zonal perturbations ( n x = 0) (panel (b)) and for non-zonal perturbations 
with n x = 1 (panel (c)) and n x = 2 (panel (d)). The non-zonal n x = 1 eigenfunction is 
the most unstable mode while the jet is stable to zonal jet (n x s= 0) perturbations, (e): 
Demonstration of the convergence of the power method iteration towards the maximum growth 
rate for perturbations with n x — 1,2. Note that the coefficient of linear damping of the mean 
flow (whether zonal or non-zonal) is ten times smaller than that of the perturbation field. The 
forcing is NIF (detailed parameter specification can be found in Fig. 6.2). 



Figure 6.10: The growth rate of the maximally growing S3T eigenfunction of the 1-jet 
equilibrium shown in Fig. 6.9a as a function of the zonal wavenumber n x of the perturbation 
(panel (a)) and similarly for the 2-jet equilibrium at e/e c ,z = 800 (panel (b)). Parameters as in 
Fig. 6.2. 


117 






















as function of the n x is shown in Fig. 6.13. While the 3-jet is stable to jet 
perturbations it is unstable to non-zonal perturbations with n x < 7. (Note that 
the stochastic forcing has power only at (7,0)). Substantial growth occurs for 
n x = 1,3,6 and maximal growth for n x = 6. While there is instability at both 
n x — 3 and 6 at finite amplitude all instabilities are attracted to a n x = 1 state 
with maximal power at n y — 4 with approximately the spectral structure of 
the n x — 1 instability. This phenomenon of saturation of the instabilities to 
the structure of an instability of smaller growth rate requires further study; it 
appears that there is at play an Eckhaus instability for non-zonal states. We 
witnessed similar behavior also in cases in which the homogeneous state was more 
unstable to non-zonal perturbations but the instabilities saturated into zonal jets. 
Thus the transition from a zonal to a non-zonal flow in the S3T simulation of 
section 4.3 is caused by the secondary S3T instabilities of the finite amplitude 
3-jet saddle equilibrium that emerges from the homogeneous background. This 
complex of phenomena predicted by S3T is remarkably reflected, with the same 
time sequence, in corresponding NL simulations (see Fig. 4.4 in section 4.3). 


118 


(a) 


* e (y) 


(b) 


n x = 0 , s = —0.0097 


0.02 

o.oi 


-o.oi 

- 0.02 


(d) 


n x = 2 , s = 0.006 + 0.358i 


(g) 



( e ) 


n x = 3 , s = 0.025 + 0.605i 


n x = 5 , s = 0.015 + 0.547i 

’#*###***#< 

• • • • 

If IfIff If 



.. 

[(((((((((((I 

•0000S000001 



Figure 6.11: Stability of the zonal jet S3T equilibrium for the case presented in Fig. 4.8. (a): 
Contour plot of the equilibrium mean flow streamfunction, T e (y), with the zonally averaged 
zonal velocity, U e (y ), (thick black line), (b-h): The most unstable mean flow perturbation 
streamfunction, for perturbations with n x = 0,1,..., 6 together with their corresponding 
eigenvalues, s. While this zonal jet S3T equilibrium is stable to zonal jet perturbations (n x = 0) 
it is found to be unstable to non-zonal perturbations with maximum instability occurring for 
n x — 1,3,6. The maximal growing eigenfunction has the (1,4) structure that eventually emerges 
both in the NL and in the generalized S3T simulation (cf. Fig. 4.4 and Fig. 4.7 respectively), 
(i): The evolution of the growth, A, for n x — 1,2 perturbations showing the convergence to 
the corresponding growth rates, s r (dash-dotted line). Parameters: = 4.9, linear damping 
coefficient r — 0.01, stochastic forcing with single harmonics with wavenumber (7,0) and energy 
injection rate £ = 4 x 10 -5 . 
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Figure 6.12: (a) The energy spectrum of the jet equilibrium shown in Fig. 6.11a as a function 

of the meridional wavenumber n y . (b-h) The energy spectrum of the mean flow perturbations 
with n x = 0,1,..., 6 as a function of the meridional wavenumber n y . The n x — 0 eigenfunction 
is a q y = 0 Bloch state, while eigenfunctions n x — 1,...,6 are q y = 1 Bloch states. Other 
parameters as in Fig. 6.11. 



Tlx 

Figure 6.13: The growth rate of the maximally growing S3T eigenfunction of the 3-jet 
equilibrium shown in Fig. 6.11a as a function of the zonal wavenumber n x of the perturbation. 
Other parameters as in Fig. 6.11. 
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6.3 Stability of finite amplitude non-zonal traveling wave S3T 

EQUILIBRIA 

Up to now we have investigated the stability of zonal equilibria to zonal and 
non-zonal perturbations. We now show how to calculate the stability of non-zonal 
wave S3T states. For the calculation we adopt the method for calculating the 
sensitivity to initial conditions of a trajectory (Z(x, t), C(x a , x&, t)) of the S3T 
system. If the state trajectory is perturbed by (SZ, SC) then its linear stability 
(i.e. its Lyapunov exponent) is obtained from a simultaneous integration of the 
system: 


dtZ + J (T, Z + (3 • x) = K(C) - Z , 
dtCab = [A(U) + M U)] Cab + eQ ab , 
dt SZ = ^(U) SZ + K(8C) , 
dt 8C ab = [A(U) + A(U)] 8C ab + (< 5A a + SAt) C ab . 


(6.5a) 

(6.5b) 

(6.5c) 

(6.5d) 


While with this method we determine unequivocally the largest Lyapunov exponent 
of a perturbation trajectory, we can also estimate the stability of a traveling 
wave state of the S3T system if the S3T state persists in that state despite being 
unstable. This often occurs because the unstable S3T states have many stable 
directions (they are saddles) and because the S3T equations are noiseless the 
unstable directions take a long time before they obtain appreciable magnitude. 
Under these conditions integration of the tangent linear system (6.5c)-(6.5d) 
produces good estimates of the growth rate of the unstable traveling wave solutions 
of (6.5a)-(6.5b). 

With this method we consider the stability of the traveling wave S3T state 
shown in Fig. 6.14. This finite traveling wave state emerged as a quasi-equilibrium 
of the finite amplitude equilibration of the most unstable non-zonal eigenfunction 
of the homogeneous equilibrium. We show that this non-zonal state is unstable 
to zonal perturbations {n x = 0) which saturate into a predominantly zonal state. 
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Figure 6.14: Stability of a non-zonal traveling wave quasi-equilibrium S3T state, (a): A 
snapshot of a finite amplitude state streamfunction, T, that consists of a wavenumber (1,4) 
wave traveling westwards. This state after long time (around 3000 time units) transitions to a 
wavenumber (0, 2) nearly zonal jet structure, (b): The mean flow streamfunction, (5T, of the first 
Lyapunov vector of the state shown in (a) together with its corresponding Lyapunov exponent, 
A. The Lyapunov vector of the doubly periodic quasi-equilibrium is neither monochromatic in x 
nor y. However, its maximum power is concentrated at (n x ,n y ) = (0,2). (c) Demonstration of 
the convergence towards the maximum growth rate, A = 0.0033. Parameters: [3 = 12, r = 0.01, 
v m = 0.005, IRFn forcing with kf = 6, 8kf = 1 and s — 3£ Cj nz = 5.24 x 10 -4 . Stability 
calculations were performed at resolution N x = N y = 32 which results in a state variable of (6.5) 
with dimension 2 x 10 6 . 
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Conclusions 


In this thesis we have presented a theory for the formation and maintenance of 
eddy-driven jets in planetary turbulence. The theory presented, called Stochastic 
Structural Stability Theory (S3T), is a non-equilibrium statistical theory which has 
its basis in wave-mean flow interaction theories. It studies the closed dynamics of 
the first two cumulants of the full statistical state dynamics of the flow, and neglects 
or parametrizes third and higher-order cumulants. Neglect or parametrization 
of third and higher order cumulants is equivalent to neglect or parametrization 
of the eddy-eddy interactions in the equations of motions. Thus, in the specific 
closure only the non-local in wavenumber space interaction between large-scale 
flows and the smaller scale eddies is allowed while local wavenumber interactions 
among eddies are not allowed. We study within this statistical closure large- 
scale structure formation and maintenance in stochastically forced-dissipative 
barotropic turbulence on a /3-plane. 

We have performed an extensive comparison of jet formation as predicted by 
the S3T dynamics and as predicted by direct numerical simulations of the fully 
nonlinear equations. The emergence and equilibration of zonal jets in the S3T 
dynamics and, moreover, the remarkable agreement of their predicted structure 
with the jets observed in direct numerical simulations, provides a constructive 
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proof that turbulent cascades are not required for the formation of zonal jets 
in /3-plane turbulence. Emergence and equilibration of zonal jets occurs due to 
cooperative quasi-linear mean flow-eddy interactions that is captured by S3T. 

The S3T dynamics is autonomous and deterministic and therefore may have 
fixed point solutions. These solutions are statistical equilibria of the turbulent 
flow that describe both the large-scale structure mean flow (1st cumulant) as well 
as the second-order eddy statistics (2nd cumulant). S3T allows the study of the 
stability of such equilibrium states. Instability of a statistical equilibrium signifies 
transition of the turbulent regime to a different attractor (i.e. to a different climate 
state). Such an equilibrium state of the statistics is the homogeneous turbulent 
state with no mean flow. We have demonstrated that the homogeneous state 
becomes unstable at analytically predicted critical parameter values and the 
flow undergoes a bifurcation becoming inhomogeneous with the emergence of 
large-scale zonal and/or non-zonal flows. The mechanisms by which the turbulent 
Reynolds stresses organize to reinforce infinitesimal mean flow inhomogeneities, 
thus leading to this statistical instability, are extensively studied for various 
regimes of parameter values (planetary vorticity gradient, dissipation rate and 
turbulent energy injection rate). It is shown that for small and modest values 
of the planetary vorticity gradient, /?, the upgradient fluxes responsible for the 
large-scale structure formation instability are induced by the Orr mechanism, 
while for large /3 they are induced by resonant wave triads. The dependance of 
the instability on the spectrum of the stochastic excitation was also studied. 

The relation between the formation of large-scale structure through modula- 
tional instability and the S3T instability of the homogeneous turbulent state was 
also investigated. We have demonstrated the formal equivalence between the 4MT 
system, that approximates well the modulational instability of coherent Rossby 
waves, and the S3T instability of a homogeneous turbulent state. However, we 
have demonstrated that the 4MT dynamical framework is inadequate for capturing 
the finite amplitude equilibration of the instabilities. 

We also presented methods for studying the stability of inhomogeneous turbulent 
equilibria. We demonstrated that the phenomenon of jet merging is properly 
understood as an S3T instability of a finite-amplitude mean flow state and 
therefore is properly understood in the framework of statistical state dynamics. 
Also, we have shown that the transition from zonal to non-zonal turbulent states 
is also predicted by S3T stability analysis. 

The S3T closure produces an analytical, predictive and quantitative theory for 
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turbulence that proceeds directly from the equations of motion. It provides a way 
of determining turbulent statistical equilibria (climate states of our model) and 
moreover determining their stability. 
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A 

Construction of the stochastic forcing and 
demonstration that the energy injection rate 
induced by a temporally delta-correlated 
forcing is independent of the state of the 

system 


The stochastic equations (i.e. eqs. (2.1) or (2.6)) are nonlinear differential equations 
with additive stochastic excitation of the generic form dt4> — £(</>) +A f(4>, (j)) + 
with C a linear operator and J\f a nonlinear operator. The stochastic function 
£(x, t) is taken to be a Gaussian process with zero mean, i.e., (£(x, t)) = 0, and 
a homogeneous function of both space and time. It is further assumed to be 
delta-correlated in time but spatially correlated with spatio-temporal covariance 
(£(x a , £)£(x&, £')) = Q(x a , x^)5(t — t f ) (the averaging operator (•) denotes an 
average over forcing realizations). The delta correlation in time is a very crucial 
assumption because it implies that the energy injection rate is constant and 
independent of the state the system and depends only on the amplitude factor e 
when Q is appropriately normalized. This allows us to know the energy injection 
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rate once we have specified the forcing without any reference to the flow that 
develops. If the stochastic excitation were not delta correlated, i.e., if for example 
it were a red-noise process, the energy injection rate would depend on the state 
of the system. 


A.l Construction of the stochastic forcing 


We show first that the spatio-temporal homogeneity of the random field £ implies 
that (£(x a , £)£(x&, t')) = F(x a — x&,£ — £'), which in turn implies that the spatial 
covariance of the forcing, Q, is a function of the difference coordinate x a — x&, 
i.e., Q(x a , x b ) = Q(x a - x b ). 

Since £ is chosen to be homogeneous its statistical properties are invariant 
under translations, that is (£(x a , £)£(x&, £')) = (£(x a + a,t)£(xb + a,t')) for any 
a,. By expanding £ in Fourier series, 

COM) = J ^^l(M)e lk ' x , (A.l) 

with £(—k,t) = £(k, £)* in order that £ be real valued, we have that: 


(CO^MOM')} 



d 2 k d 2 k' 
( 27 r ) 2 ( 27 t ) 2 




gik-Xagik'-x^ 


(A.2a) 


while 


(C(x 0 + 0C(x& + a, t 1 )) = Jj 


d 2 k d 2 k ; 
( 27 t ) 2 ( 27 t ) 2 


I(k, t)i( k', t')) g ik ' x <i e ik '- x 6 gi(k+k')-a 

(A.2b) 


For (A. 2a) and (A. 2b) to be equal we must have e 1 ( k + k, )-“ = 1 for every a, which 
requires that k 7 = —k and which in turn implies that the wavenumber covariance 
of the field is of the form: ^£(k, £)£(k', t')^ — ( 27 r ) 2 ^£(k, £)£( —k, t')^ 5(k + k'). 
Consequently, 

<C(x a , *)C(*6, 0 ) = / (|( k , i)|(-k, 0) e ik -( x “-^) , (A.3) 

which is a function only of the difference x a — x&. Similarly, we can show that 
temporal homogeneity implies that the ^£(k, £)£(—k, = /(k,£ — £'). 

Note if the random field is homogeneous only in the zonal x direction, then 
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the above argument would require that = \ for every a x , which 

would imply k x = —k! x and as a result the spatial covariance of a random field 
homogeneous in x should be a function of x a — x b , y a and y b . 

To construct the delta-correlated Gaussian stochastic excitation we choose the 
Fourier amplitude of (A.l) at each wavenumber k and at each instant to be: 

l(M) = w(k)r] Kt , (A.4) 


with w( k) the amplitude at wavenumber k, satisfying w(— k) = w( k)*, and r/k ,t 

a complex valued Gaussian random variable of k and t satisfying 

with zero mean and delta-correlated covariance both in wavenumber and time: 

t') ~ — t')5(k — k'). Equation (A.4) implies that (£(x,£)) = 0 and, 

moreover (A.3) can be written as, 


(£(Xa,*)£(Xb>0) = JJ 


(2?0 2 (LI ) 2 %',*') e lk ' Xa e lk '- X6 

= // » (k M k 'r (%,«4,,,,) ■=» 


= S(t-t') j , 


(27r 


(A.5) 


(in the second equality we changed the integrating variables k' —>► —k'). Therefore 
the spatial structure of the covariance, Q, is given by: 


Q(x a - x b ) = 


/ 


d 2 k |^(k)| 2 
(27 t) 2 (27t) 2 


Q ik-(x a —Xfr) 


/ 


d 2 k 

W ) 2 


Q(k) e 


,ik-(x a -x b ) 


(A.6) 


The spatial covariance Q turns out to be an even function of its argument, since 
| w HOI = |u;(k)|. This is physically expected due to the homogeneity of £, which 
implies the exchange symmetry (£(x a , £)£(x&, t')) = (£(x&, t)£(x a , £')). Moreover, 
equation (A.6) shows that the Fourier transform of Q besides being real is also 
non-negative, as is demanded by Wiener-Khinchin theorem. The statement that 
homogeneous covariances have real and positive Fourier transforms is often also 
referred to as Bochner’s theorem. 

That the covariance (£(x a , £)£(x&, t r )) turns out to be a function only of the 
difference coordinates x a — X5 may be alternatively seen as a result of the fact 
that each of the Fourier components of the stochastic forcing corresponds to a 
different random variable which is uncorrelated to the others, except for the pairs 
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k and — k. 

We construct a discrete representation of £ in time as follows. Discretize time 
with time steps h. The complex random variable at time t = jh, j = 0,1,..., 
is taken to be 

Vk,t = ^k,t) •> (A.7) 

with X and Y real valued random numbers taken from a Gaussian number 
generator with zero mean and unit variance; we further set X_\^ t — X^ t and 
= ~Yk t . The Vh in the denominator ensures that (A.7) is delta-correlated in 
time satisfying as h -A 0 the discrete expression of the delta function requirement, 
Ho (rjk dt' = 1. Note that because X and Y are normally distributed, the 
probability density function of 77 k,t, T* (%,*)> is only a function of the amplitude of 


P(Vk,t) = P(X Kt ,Y Kt ) 

= P(X Kt )P(Y Kt ) 

1 _I y 2 1 _lv2 
— _p 2 k,i _^ 2 J k ,t 

= — e ~^M 2 , (A.8) 

2n 

which is only a function of the absolute value of the 77 k, t- 

A. 2 Proof of the relation (£(x a ,£)C'(x 6 , t) + C'(x a ,£)f(x & ,i)) = y / eQ (* a - 

In this section we calculate the term (£(x a , £)£'(xb, £) + £'(x a , £)£(x&, £)) i n (2.11), 
where £' satisfies the NL eddy equation (2.4b). 

We should note that because the stochastic forcing is additive (it does not 
depend on the state of the system) there is no need to distinguish between 
the mathematically convenient Ito interpretation of the stochastic differential 
equations, in which the state at t is uncorrelated with the stochastic forcing at 
the same time £, i.e., (^(x, £)£(x, t)) — 0, and the physically relevant Stratonovich 
interpretation in which state and forcing are correlated at the same time. Both 
interpretations lead to the same results (0ksendal (2000, p. 35)). Throughout 
this thesis we have adopted the Stratonovich interpretation. 

To calculate (£(x a , t)(’ / (x 5 , t) + £'(x a , £)£(x&, t)) we write the solution of (2.1) 
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in integral form as 


COM) = C(x,0) - J J (V'(x,s),C(x,s) + /3 • x) - rC(x,s) ds + Vi j £(x,s)ds. 

(A.9) 

from which it follows that the eddy vorticity is given as: 

C'( x ,*) = C'( x ,0) - |j (^(x,s),C(x,s)+^ x) 

J (^(x,s),C(x,s) + /3 -x) - rC'(x, s) j ds + V^j^ £(x,s)ds . 

(A-10) 


- T 


We want to calculate (£(x a , t)(’ / (x 5 , £)). From (A.10) we have: 


(^(x a ,t)C'(x 6 ,t)) = 

= (£( x a,0C'( x b,0)) - J q ({£( x a,0 J (V’( x 6,s),C( x 6,s) +/3 • x) 

- £(xa,t)T J (tp(xb,s),C(x b ,s) + /3 -x) - r£(x a ,f)C'(x fc ,s)}^ ds 

+ V? [ (£( x a,i)£( x fe, s )) ds . (A.11) 

Jo 


The initial state of the system is clearly uncorrelated with £ at time £, (£(x a , £)£'(xfc,0)) 
0. The first integral on the r.h.s. does not contribute to the correlation since 
i) there is no correlation between stochastic forcing at time £, £(x, £), and state 
of the system at any time t' < t and ii) also because the integral is unchanged 
if calculated over the interval s G [0, t) instead over s G [0, t] since at time £, 
(£(x a , t)(’ / (x 5 , £)) < oo. The excluded point s = £ is only a point of measure zero. 
Therefore, the only non-zero contribution to (£(x a , £)£'(xb, £)) is the last integral, 


<^( x a ,t)C'( x 6,i)) = / (£(x a ,i)C(x&,s)) ds 

Jo 

= Vs g(x a - Xfe) (5(t - s) ds 
JO 

= “y C( x a - x 0 , 


(A.12) 


where the integration of the delta function gives 1/2. Similarly, (£'(x a , i)£(x&, t)) = 
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(V^/2)Q(x a - x 6 ), leading to 

(f(x a ,i)C'(x 6 ,£) + C'(x a ,£)f(x & ,£)) = VzQ(*a ~ X 6 ) . (A.13) 

A. 3 Energy injection rate in the NL and QL systems by the stochas¬ 
tic EXCITATION 

The energy of the flow is defined as E = / d 2 x ^|u| 2 and the time rate of change 
of E, assuming periodic boundary conditions, is given after an integration by 
parts as 

if = T / ^ = ~Jd 2 xrpd t C . (A-14) 

We calculate the contribution of each term of dt( from (2.1). After integration 
by parts, the Jacobian term gives 

J d 2 x 'ip J(^, ( + /3 • x) = 0 , (A.15) 

as expected because the Jacobian terms merely redistributes the energy among 
the flow scales. The rate of energy due to the dissipation is 

r J d 2 x = — r J d 2 x |V?/;| 2 = —2 rE < 0 . (A.16) 

The ensemble average energy injection rate from the stochastic forcing is 

-Ve J d 2 x (V>0 • (A.17) 

Since £ is a stochastic variable we seek to determine the ensemble average energy 
injection rate over all forcing realizations. To obtain this estimate we proceed 
as in (A.9). As argued above, the first integral on the r.h.s. of (A.9) does not 
contribute at all to the correlation (^(x, £)£(x, t)) and the only contribution comes 
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from the last term in (A.9) so that, 


(£(x, t)ip(x,t)) = (£(x a ,i)^(x 6 ,i)) 


X a =Xfe 


= ^a l Q{*a-Xb) 


Xa=*b 


= Vl f d2k 0( k ) „ik-(xg—: 


[ a «- VW ik-(x„-x 6 ) 

J (2t r) 2 -fc 2 

d 2 k Q(k) 


x ri =x fe 




(27t) 2 2fc 2 

From (A. 17) we have that the energy injection rate is 


- v 7(/d I x,Kx, 1 )«x 1 t))= £ ^?ffi) (/d 2 x) 

Therefore, if the covariance spectrum is normalized according to 

[ d 2 k Q(k) _ 

J (27r) 2 2fc 2 


(A-18) 


(A-19) 


(A.20) 


then the energy injection rate per unit area is e. 

From (A. 15), (A. 16), (A. 19) and (A.20) we conclude that the domain and 
ensemble averaged total energy of the flow, (. E) = (/ d 2 x) -1 F), satisfies: 

^=£-2r (E) , (A.21) 

and this implies the total energy the flow (the sum of the eddy and mean flow 
energy) will always approach with time the constant value: 

(^)oo = ^ ^ 
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B 

General properties of S3T equilibrium 

solutions 


Assume that (Z e ,C e ) is an equilibrium solution of the S3T statistical state 
dynamics governed by (2.13) satisfying: 


J (^ e , Z e + (3 • x) - K(C e ) + Z e = 0 , (B.la) 


A(u e ) + A(u e ) 


Cab + £ Qab ~ 0 , 


(B.lb) 


with U e = z x VA ~ x Z e . 

We show first that when A e = A(XJ e ) is stable then there exists a unique 
Hermitian and positive definite solution C e of (B.l). 

When the operator (, A e a + AD is invertible then (B.lb) has a unique solution. 
The spectrum of (, A e a + AD is Hi + /i*, i,j = 1, 2,... where Hi are fh e eigenvalues 
of A e . If A e is stable then R e(Hi) < 0, which implies that R e(ni + Hj) < 0 for 
every z, j and hence (A e a + AD is invertible. So a unique C e that solves (B.l) 
exists. We now show that it is a physically realizable solution, i.e., it is Hermitian 
and positive definite, by giving the explicit expression of the solution. We verify 
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(B.2) 


that the solution of (B.l) is 


C e (x a ,x & ) 


e lim 

t—>•+oo 


/ 


AH. 


v*»- 


! g(x a - x fe ) ds 


since 


(a® + .4g) = £ t Hm^ ^ (a® + Afj e A ° s e^ s Q a& ds 


= £ lim / * - 

t —^-|-oo 


= £ lim e A * s e A t s Q ab 

t —^-|-oo L J s — 0 

= -£ Qab + £ lim e A °- t Q ab 

t^+OO 


it 


£ Qab i 


(B.3) 


as lim^ +00 e^ t e' A b t Q ab vanishes when A e is stable. 

From (B.2) it is clear because Q(x) = Q(— x) that C is Hermitian, i.e., it 
satisfies the exchange symmetry C a b = Ci a . It is only left to prove that C e is 
positive definite. A real function F(x a ,x^) is positive definite if and only if for 
every complex function /, 

JJ dx a dxb F(x a ,x b ) /(xo)/(x fe )* >0. (B.4) 

Consider a positive definite homogeneous covariance, like the forcing covariance 
Q. Then for all functions (or even distributions) / we must have: 


JJ dx a dxft Q(x a - x 6 ) /(x a )/(x 6 )* = 


= / JJA JJ dx a dx fe Q(k)e lk ' (x “ Xfc) /(x a )/(x 6 ) 

= /(!l^« k) - /<^)) (/ 


dx fe e ik ^/(x 6 ) 


-I 


d 2 k 

(2vr) 2 


Q(k) /(k) >0, 


(B.5) 


This implies that Q(k) > 0 given that this integral must be positive for all 
functions /. (This is Bochner’s theorem.) We now show that C e is also positive 
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definite. Indeed from (B.2) we have for all functions /: 


JJ dx a dx fc C e (x a ,x b ) f(x a )f(x b )* = 


= E 


// dx “ dxt s: 


d s 


e -4a S e A e b S Q^ Xa _ f ^x a )y*(x^)* 


= e l iV \f 7^rp IJik > /dx (e^* V k x )/(x) 


e^ s e^ s Q(k)e ik ' (Xa “ Xi>) 
2 


/(x a )/(x 6 )* 


> 0 . 


(B-6) 


Note that, 

C(x a , x fe , i) = e- 4 ^ t e* 4 6 *C(x a , x 6 , 0 )+£ [ e A “ s e A b s Q(x a - x b ) ds , (B.7) 

Jo 

solves the time dependent S3T covariance equation (2.13b), 

dfCab = {A e a + A^j Cal) + £ Qa6 • (B.8) 

To show that make the change of variable s t — s in (B.7) and write the 
covariance in the form 

C(x a , x b , t ) = e A “ t e A b t C{x a , x b , 0 ) + £ f e - A “ ^~ s ^e A b (t_s) Q(x a - x b ) ds . 

Jo 

(B.9) 


which satisfies (B.8). From (B.7) we see that if C(x a ,Xb, 0) is positive definite, 
then C(x a ,Xft,t) is positive definite at all times £, irrespectively of whether A e 
is stable or unstable and when A e is hydrodynamically unstable (7(x a ,x^,t) 
diverges as t oo. This implies that no S3T equilibrium exists when U e is 
hydrodynamically unstable. 

It should be noted that the time independent Lyapunov equation (B.lb) has 
always formally a solution unless A e has a neutral eigenvalue, i.e. unless there 
is an eigenvalue with Re(/i^) = 0. However, the formal solution of (B.lb) for 
A e unstable is not positive definite and consequently is not physically realizable 
given that all physically realizable states produce positive definite covariances (or 
from the previous argument this non-positive steady state covariance solution is 
unreachable from any physical initial covariance which generates through (B.9) 
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only positive definite covariances). In order to see that consider the eigenrelation 
of A e and its adjoint A e V 

A u n = ii n u n , A ^ v n = ii n v n 5 72 = 1, 2,... . (B.10) 


The eigenvalues of A e ^ are the complex conjugate of the eigenvalues of A e and in 
general the eigenfuctions u n of A e and eigenfuctions v n of A e ^ are not the same. 1 2 
Moreover, they satisfy the following relations: 


\Um,v n ) = 6 mn (u m ,v m ) (orthogonality), (B.lla) 

y, °^ti( = S(x a — Xft) (completeness) , (B.llb) 

n y^n-> u nj 

where (•, •) denotes the inner product between two functions g and /i, taken here 

as ( g,Jij =/d 2 xjf(x)*ft(x). 

Assume that A e has an eigenfunction, say u±, with eigenvalue with Re(/ii) > 0 
with associated adjoint eigenfunction v\. We will show that the solution C e 
of (B.lb) is not positive definite. Consider the integral 


J d 2 x a J d 2 x 6 A e a C e {x a ,x b )\ vi(x a )v!(x b ) 

= f d 2 x a f d 2 x 6 f dV A e a y 

j n [v n ,U n j 


MXaK(x') c e (x / ; ^ ^( Xa ) Ul(x6 ) 


= J d 2 x fe J d 2 x' y Mn 


(^Vi , Uji^j 


C e (x^X 6 )<(x>i(x 6 ) 


Mi J d 2 x a J d 2 x b C e (x a ,x 6 )^(x a )ni(x 6 ) , 


(B.12) 


and similarly, 


/ d2x “/ 


d 2 x fe ^C e (x a ,x 6 ) ni(x a )ni(x fe ) = 


= Mi J d 2 x a J d 2 x b C e (x a , Xft) n^(x a )ni(xfc) . 


(B.13) 


1 The adjoint of operator .4 e is defined as the operator that satisfies 

f d 2 x g(x)* [*4 e fi(x)] = f d 2 x [A et g(x)j* fi(x), for every functions g, h. 

2 Eigenfunctions u n and v n coincide only when the two operators commute, i.e., A e A eif = 
A ef A e , in which case the operator A e is called normal. 
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Using (B.12), (B.13) and (B.lb) we have that 



(B-14) 


hence, because Re(/ii) > 0, we conclude that 



(B.15) 


which shows that C e is not positive definite and hence not realizable when A e is 
unstable. 

Further, in the absence of forcing, i.e., Q = 0, the equilibrium solution (Z e , C e = 
0) is hydrodynamically stable if and only if it is S3T stable. To see that consider the 
S3T perturbation equations (2.25) about (Z e ,C e = 0) which for this equilibrium 
simplify to: 


dt SZ = A e 6Z + K(8C) 


(B.16a) 

(B.16b) 



From the above equations we infer the equivalence between the stability of A e and 


the S3T stability of (Z e ,(7 e ), which is governed by the upper-diagonal operator 



(B.17) 
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c 


Numerical integration of the NL, QL and 

S3T systems 


In this appendix we present details of the numerical method used to integrate 
the NL system (2.1), the QL system (2.6) and (2.17) and the S3T system (2.13) 
and (2.21). The flow domain is taken to be a rectangle of size L x x L y , with 
periodic boundary conditions on its boundary and discretized with N x x N y points, 
N x in the zonal and N y in the meridional direction. The allowed wavenumbers in 
this domain, for N x , N y even, are: 



(C.l) 


'3 


and moreover, the continuous Fourier transforms become discrete Fourier trans¬ 
forms: 



(C.2) 
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with the continuous and the discrete Fourier amplitudes (j>{k x ) and (pf ix related by 


N x / 2-1 

4>(k x ) = 2ir ^2 4>k' x S(k x - k' x ) . 


k' x =-N x / 2 


(C.3) 


To simplify notation we denote: 

Nj/ 2—1 

E > j = x,y> (c.4) 

^=i 


for any variable L 


JVj/2-l 


E = E « E = 


tj=-Nj/2 

kj^O 


C.l Numerical integration of NL and QL system 

Both NL and QL systems are a stochastic partial differential equations of the 
form: 

t) = C (</>(x, t)) + J\f (V>(x, t),<p(x, t)) +y/e^{x,t) , (C.5) 

with C a linear operator and J\f a bilinear operator. We use pseudospectral 
methods to evaluate £(</>) and Nicj), (/))• Time is discretized in equal steps, /i, and 
the deterministic part of (C.5) is advanced in time by h using a fourth-order 
Runge-Kutta (RK4) scheme. After the completion of each RK4 step, the stochastic 
excitation term, hy/£%(x, nh), is added to the state of the system (Godunov, 1959). 

C.2 Numerical integration of the S3T system 

In this section we discuss the numerical integration of the S3T system. While 
the state of the NL and QL has dimension N x N y , the S3T system has the far 
larger dimension N x N y + N^Ny. This makes the numerical integration of the 
S3T system extremely costly and special techniques need to be developed in order 
to obtain resolved simulations. 


C.2.1 Numerical integration zonal S3T system 

The S3T system which is computationally more easier to integrate is the zonal 
mean S3T system (2.21) because of the homogeneity in the x direction. In that 
case the dimension of the mean flow is N y and the dimension of the covariance is, 
as we will show, at most N x Ny. 
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We first use the homogeneity in x to split the covariance equation (2.21b) into a 
system of decoupled equations for the zonal Fourier components Ck x of (7, defined 
as 


C( x a, *b,t) = E Ck x (: y a , Vb, t) e lkx ( Xa Xh) , (C.6) 

kx 

while the spatial forcing covariance is similarly expanded as 

Q(x a - x b ) = (3/a - Vb) e ikx ( Xa ~ Xb ^ . (C.7) 

kx 

Note that the zonal Fourier component Qk x is related to the forcing spectrum 
Q{ k) = / d 2 (x a - x b ) Q(x a - x b ) e - ik -(x a -x l ,) ) w ith k = (. k x ,k y ), by 

Oh(ya-y6) = Ee( k ) eifc!/(2/a “ 2/b) • ( c - 8 ) 

ky 

Note also that if we take the x Fourier transform of 

(C.9) 

kx 

the eddy vorticity covariance, being homogeneous in x, can be written as 


C(x a ,x b ,t) = ( C\x a ,t)C'(x b ,t )) 

= E (c k x (ya,t)( k , x (y b ,t)) e ikxXa e ik 'x Xb 

kx,k' x 

= E ( Ck x (ya,t)Ck x (y b ,tr) e ikx ^~ Xb ^, (C.io) 

kx 

which implies that 

Ck x (y a ,y b ,t) = (Ck x (y a ,t)Ck x (yb,t)*) ■ (C.ii) 
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We can write: 


A z ,a{U)C(x a ,yL b ,t) = 

= E [ _ ik x U * - ik * (p - d y a ya U °) A a 1 - l] C kx (y a , y b , t ) e ik A x a-x h ) ? 
kx 

= E A z , kx , a (U) C kx (y a , yb, t) e ik A x «~^ , 


where 

A,fc B (^) = - ifc* (/? - ^z* 7 ) A” 1 - 1 , 

and similarly 

A z , b (U) C(x a , X6, t) = E A z ,k x ,b(U)*C kx (y a , y b , t) e ifc « (*““**> . 

kx 

Combining (C.6), (C.7), (C.12) and (C.14) we write equation (2.21b) as: 


(C.12) 


(C.13) 


(C.14) 


J2\dtC kx {y a ,yb,t) - 


Az : kx,a(U) + *A z ,k x ,b(U) Ck x (y a ,yb,t) 

- £Qk x (y a -yb)}e l 


ik x (x a -x b ) _ q 

(C.15) 


which implies that the covariance equation decouples to the N x equations 


dtC kx (y a ,yb,t ) = 

= A Zt k x ,a(U) A Zt k x ,b{U) C kx {y a , y b ,t) + £ Qk x (ya — yb) ■ 


(C.16) 


Covariance C is real valued, implying that C-k x (y a ,ybit) — Cfc x (2/a,2/&,£), and 
therefore we only need to solve for k x > 0. Note also that the exchange symmetry 
C(x a ,x b ,t) = C(x fc ,x a ,t) implies C- kx (y a ,y b ,t) = C kx (y a ,y b ,t ) and therefore 


C kx {y a ,Vb,t) = C kx (y b ,y a ,t)* . (C.17) 


Moreover, from (C.16) we can see that we can integrate only the components k x 
for which Qk x has power (as from uniqueness if Qk x = 0 and initially Ck x = 0 
then Ck x — 0 for all times). By removing the wavenumber components for which 
Qk x — 0 we may remove also self-sustained turbulent states, i.e., the possibility 
that a non-zero covariance component is maintained in the absence of forcing 
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(which can happen when U produces a neutrally stable A Zi k x (U), when U is 
time independent, or a A z ,k x {U (t)) with zero top Lyapunov exponent when U is 
time-dependent). These self-sustained states are singular in barotropic turbulence 
because barotropic turbulence does not self-sustain and including them may 
result in singular states that disappear with the introduction of even the slightest 
forcing (see for example the states in Marston et al. (2008)). However, there 
are particularly interesting self-sustained states in 3D turbulence (cf. Farrell and 
Ioannou (2012); Constantinou et al. (2014b)). We do not include in the S3T 
calculations wavenumber components that are not externally forced. In this way 
we consider in (C.16) only the k x = 1, 2,..., N k for which Q kx is non-zero. 


We now formulate (C.16) on the discretized domain. 

The Fourier coefficients (k x (y,t) on the discretized domain become a vec¬ 


tor Ck x (i) with elements [C k x (t)\j ~ Ck x (yj-> t), while the Fourier coefficients 
Ck x (y a ,yb,t) and Qk x (y a - yb) become matrices C kx (t) and Q kx with elements 
= C , k Jy l ,y J ,t) and [Q kx \ij = Qk x {Vi ~ Vj) respectively. Also ,the 
mean flow U(y,t) becomes a vector U{t ) with elements [U(t)\j = U(yj,t). The 
symmetry property (C.17) requires that the C k and Q k matrices are hermitian, 
* - Cfc. Relation (C.ll) implies that [QJ^- = ([C/Ji or simply 


i.e., C k = 

4 


= \Ck x Ck x )i demonstrating that the C kx matrices are also positive definite. 
We express the S3T covariance equations (C.16) in terms of matrices Ck x and 
Qk x Term A z ,k x ,a( U ) + Az, kx ,b{U)* C kx (y a , yb, t) becomes in matrix notation 


[A z ^k x ,a(U) + A M (U) ] Ck x (?/a, 2/6, t) 


*k x (U)C kx (t) + C kx (t)A kx (U)l , 

J ab 

(C.18) 


where A kx (U) is the time-depended matrix 


A kx (U) = —ik x \J 


ik x 


Pi ~ (U w ) 



(C.19) 


with I the identity, A kx the matrix approximation of d yy — k A^ 1 its inverse 
and U = diag(£/), \J yy = dmg(d yy U ), where diag(#) denotes the diagonal matrix 
with diagonal elements the values of its argument at the y collocation points. We 
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thus end up with equations for ( N y x N y ) sized matrices: 


^ C-'kx (t) — A k x (U) Ck x (t) + C* B (t) A k x (E/)^ + £ , for — 1,..., JV& . 

(C.20) 

Turning to the mean flow equation (2.21a) we express the Reynolds stress 
divergence term, ^(A~ 1 d Xa +A^ l d Xh )C a b in terms of C& x . The Reynolds 

stress evaluated at y — yj is, 


X«=Xfe 




x a =x b 


2/=2/i 


_ E 9 ^ ^k Xi b) ^ a ’ ^ e ^ 


X a =Xfe 


2/=2/j 


E + c ^(*) ( ifc * 

k x 

T R e ( ( i/c * A^ 1 ) C kx (t) + C kx (; t ) (ifc* A^ 1 ) 

k x ^ 

(*)]..} • 



e i k x (x a —x b ) 1 


ab ) Xa=*b 


i/=yj 



1 


33 ) 


= E 2Re 


k x 

k x ^>0 


ik x A, 1 Cfe E 


(C.21) 


The discrete S3T system in the zonal mean-eddy decomposition takes the form: 


d. 

dt 


N k 


U(t) = ^2 2 Re vecd (ik x A k ^Ck 


k x —1 


u(t), 


^C kx (t) = A kx ( U ) C kx (t) + (t) (t7)t + e Q kx , 


(C.22a) 

for k x = 1,..., iV fc , 
(C.22b) 


where vecd(M) denotes the vector with elements the diagonal elements of matrix 
M. The homogeneity in the zonal direction has resulted in reducing the dimension 
of the S3T system from N x N y + N%Ny to N y + 2N^Ny (the factor 2 comes 
up because C^ x are complex valued). This is a tremendous reduction, i.e., for 
N x — N y — 128 and — 15 this gives around 550-fold decrease in the dimension 
of the S3T state variable. 

In previous S3T studies (for example in Farrell and Ioannou (2003, 2007); 
Constantinou et al. (2014a)) the Reynolds stress divergence term (C.21) appears 


143 



















with a factor 1/2 instead of a factor 2. The reason is that in those studies Fourier 

transforms are defined as <f>' = k x Re [4>k x e lkxX )• This results in a factor 1/2 

k x > 0 V J 

difference in the Fourier coefficients of (/ and in turn a 1/4 difference in C^ x . 

To integrate (C.22) a RK4 time stepping scheme is used. 

C.2.2 Numerical integration of the generalized S3T system 

The generalized S3T system (2.13) does not a priori posses any homogeneity in 
x, so the dimension of the S3T state is N x N y + N%N y instead of N y + 2 N^N y , 
when the flow is homogeneous in x and is represented by N\ zonal waves. 

On the discretized domain the mean flow streamfunction, 4/(x, £), is represented 
by an (N x N y )~column vector with elements: 

[*(i)]p = *(x P ,i) . (C.23) 

where index P runs through 1,2,..., N x N y covering all the (x, y ) points of the 
domain. All other mean flow fields can be expressed in terms of 4/. The eddy 
vorticity covariance C(x a ,X6,t) is represented by an ( N x N y ) x ( N x N y ) matrix, 
C, with elements: 

[C(t)]p Q = C(xp, XQ, t) , for P, Q = 1,..., N x Ny . (C.24) 

Similarly, the spatial forcing covariance is represented with matrix Q with ele¬ 
ments: 

[Q]pq = C(xp -xq) , for P,Q = l,...,N x Ny . (C.25) 

In matrix form C(x a ,X 5 ,t) = (('(x a , t)(’ / (x 5 , t)) is 

c (t) = {C(t)C(t) T ) , (C. 26 ) 

where ^ is the (N x N y )-column vector with elements 

[C(t)]p = ■ (C.27) 

We express the S3T covariance equation (2.13b) in terms of matrices C and Q. 
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For example, term ( U a d Xa + Ubd Xh )C(x. a , x&, t) becomes in matrix notation: 

[U(x a , t)d Xa + U (x b , t)d x J C(x a , x b , t) = {[U (t) D x ] C (t) + C (t) [U(i) D X ] T }^ , 

(C.28) 

with the matrix approximation of d x and U = diag (U) the diagonal matrix 
with elements the values of U (x, t) at the (x, y) collocation points. Similarly, we 
transcribe in matrix form all other terms in (2.13b), which becomes 

^C(f) = A(U) C (i) + C(t) A(U) t + £ Q, (C.29) 

with A(U) is the time-depended matrix 


A(U) = -(UD* + VD y ) + [(AU - P y \)D x + (AV + p x l)D y ] A 1 - I , (C.30) 


and A the matrix approximation of the Laplacian operator A and A -1 of its 
inverse. 

In the mean flow equation (2.13a), we express the Reynolds stress divergence 
term, 7 Z(C) in terms of the matrix C. At point xp we have: 


n(c) 

= - V- 

iz x (VaA" 1 + VfcA^ 1 ) C ab 


x=xp 

Z x ' 

= - V • 

2 

(^2/a^a 1—dxa^a 

l r 

/ i \ / i\T 

= 2 Dx 

(d,a- 1 )c + c(d,a- 1 ) 


J x«=x 6 


X=Xp 


Xa=Xfe 


X=Xp 


ab. 


+ l°y 


(d x a- 1 )c + c(d x a- 1 ) 


X a =X fe 

T 


x=xp 


- abl x a =x b 

c 


x=xp 


(C.31) 


= vecd (D y A ^ C + D y vecd — (d x A ^ 

Therefore the discrete S3T system becomes in matrix form: 

jZ(t) + J Z (t) + /3 ■ x) = R (C(f)) - Z(t) , (C.32a) 

d 


^C(t) = A(U) C(t) + C(f) A(U) t + e Q , 


(C.32b) 

where R produces the Reynolds stress divergence for covariance C and is defined 
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as 


R(M) = D x vecd 




M 


+ 


Dy vecd 


(d.a- 1 ) 


M 


(C.33) 


The Jacobian in (C.32a) is calculated pseudo-spectrally. 
RK4 time stepping scheme is used. 


To integrate (C.32) a 
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D 

Bloch’s theorem 


In this Appendix we will use Bloch’s theorem to prove (3.4) and (4.4). 


Consider the eigenvalue problem: 

*4(x) /<x(x) = C /<t(x) , (D.l) 

where A is the linear differential operator, 

-4( x ) = b(x) + b M) d Xj + 6 ifc( x )^i* fc + • • • • ( D - 2 ) 

3 3,k 

Suppose that the coefficients of A are invariant under the translation x —>> x + a, 
where a is any integer multiple of the constant vector ao- This implies that if 
f<j (x) is an eigenfunction of A so is f c r (x + a). To see that set x —x + a in (D.l) 
and use *4(x + a) = *4(x) to obtain ^l(x)/ cr (x + a) = a /^(x + a). 

Define T r the translation operator: 

T P /(x) = /(x + r) . (D.3) 
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Operators A and T a commute, i.e., AT a = T a A. This can be established by 
considering, without loss of generality, the action of T a A on an eigenfunction of A: 
[T a A(x)] /cr(x) = T a [A(x) f a (x)] =crT a / (J (x) = cr/^x + a) = ^(x + a)/ a (x + 
a) = v4(x) [T a /cr(x)] = [*4(x)T a ] / cr (x). The commutation of ^4 and T a implies 
that the eigenbasis of A and T a can be chosen to be common. 

This is not automatically achieved. For example all plane waves of the form 
e iq x are eigenfunctions of T a but not necessarily of A. From Fourier analysis 
we know that because every function can be written as a superposition of plane 
waves, the eigenfunctions of A will be of the general form f c r (x) = / d 2 q c(q) e iq ' x . 
To determine the constraints imposed by the periodicity of A on the eigenfunction 
it is only required, because A and T a commute, to render f Q r (x) an eigenfunction 
of T a . Since a is an integer multiple of ao, to obtain the most general eigenfuction 
satisfying all symmetries of A we need to render / CT (x) an eigenfunction of T ao . 
Then, since 


T ao / a (x) = Tj 


ao 


I 


d 2 q c(q) e iq,> 


/ 


= / d 2 q c(q) e iq ' 


x e iqa 0 


(D.4) 


if f (j (x) is to be an eigenfunction of T ao then e iqa ° cannot depend on q and should 
be of the form e iqa ° = e 1 ^. This holds when q = n + mpo with m any integer, 
Po the vector: 1 

Po = 27ra 0 /|a 0 | 2 , (D.5) 

that has the property that po • ao = 2n and n a constant vector satisfying 
|n| < |po|/ 2 = 7r/|ao|. This restricts the Fourier spectrum of the eigenfunction 
f (j (x) allowing only power at discrete wavenumbers. By writing c(q) in the form 


c(q) = Cm 6 (q - n - mp 0 ) , 

m 

we obtain that a general eigenfunction of T a is: 

/(x) = £ C m e i ( n + m Po)' x = e in-x ^ Cm e 


(D.6) 


impo-x 


(D.7) 


9(x) 


This implies that the general eigenfunction of T ao and A is of the form e in-x g(x), 
for any n satisfying |n| < 7r/|ao|, and for periodic functions g(x) with period ao, 


1 In solid state physics po is referred to as the fundamental vector of the reciprocal lattice. 
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i.e., g(x + ao) = g(x). This the content of Bloch’s theorem about the structure of 
the eigenfunctions of operators with periodic coefficients. 

In the special case that A is homogeneous, i.e., its coefficients do not depend 
on x, A commutes with T a for every vector a. From (D.7) we see that in this 
case g(x) must be constant and therefore the only common eigenfunction of A 
and of all the T a ’s is the single harmonic 

/cr(x) = Co e in ' x , (D.8) 

with no restriction on n since |ao| may be taken infinitely small. 

Another way at arriving at the same result is to note that eigenvalue e ina ° of 
Tao is degenerate and the degenerate eigenfunctions form a subspace G ao spanned 
by e iq ' x with q = n + mpo, m — 0, ±1,..., po • ao = 2n and |n| < 7r/|ao|, as 
previously. The common eigenfunctions of A and T ao will thus be indexed by 
n and for each n will be linear combinations of basis of G ao (as in (D.7)). To 
arrive at the homogeneous result (D.8) consider a decreasing sequence of ao, i.e., 
ao, ao/2, ao/2 2 ,..., each associated with its own degenerate subspace for an n. 
The eigenfunctions of A for that n will belong to H^i G ao / 2 s, which has the only 
element: e inx , with n unrestricted. 
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E 

Derivation of the eigenvalue relation for the 
S3T stability of a homogeneous equilibrium 


E.l Eigenvalue relation for homogeneous S3T equilibrium 

Here we derive an analytic expression satisfied by the eigenvalue a that determines 
the stability of the homogeneous equilibrium (3.1). The eigenvalue problem to be 
solved is (3.3), 


s8Z = A e 8Z + K(6C) , (E. la) 

a SC ab = [A e a + At) 6C ab + (5Aa + SA b ) C e ab . (E. lb) 

We have included in (E.l) also hyper diffusive damping of the form: —(—l) h V2h^ h C 
For h — 1 this is normal diffusion, while for h > 1 is hyperdiffusion of order 2 h. 
With this damping the operator A e becomes: 

= z • (/3 x V) A" 1 - 1 + (-l) h v 2h A h . (E.2) 
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Equations (E.l) is a linear system in 5Z and SC and can be written symbolically 
as: 1 


s 


SZ , SC = C(C e ) <SZ, 


iT 


(E.3) 


Because the only spatial dependence of the coefficients of C is through C e {x a — x b ) 
it is advantageous to consider the eigenvalue problem in terms of the relative coor¬ 
dinate x a — x b and centroid coordinate (x a + x b )/2. Since (E.la) is homogeneous 
in x and (E.l) homogeneous in (x a + x b )/2 the eigenfunctions SZ and SC will be 
single harmonics of their homogeneous coordinate (cf. Appendix D) and can be 
assumed to be of the form: 


5Z n (x) = e in x , (E.4a) 

8C n (x a ,x b ) = C(x a - x b ^j e in -(xa+x 6 )/2 . (E.4b) 

The homogeneous part of the perturbation covariance eigenfunction, C^\ is 
expanded as 

ci h) (x a - x b ) = I A!|_c( h )(k) e ik-(x a -x,) } (E .5) 

and introducing (E.5) to the perturbation covariance equation (E.lb), we obtain: 


(Aq + A b ^j SC ab — — J ^ u<2 + ^ 2 -) + i (^k + — ^k_) C^(k) 


Q ik -(Xa-Xb) 


(E.6a) 


lr To find explicitly the form of C one needs to manipulate 5A a C^ b in the following manner: 

SAa c e ab = —<5U a • V a C e ab + (A a <5U a ) • VaA" 1 ^ 

= - (z X VA'I'a) • V a Cl b + (z X V a 5Z a ) ■ V a A 
= (z x V a C e ab ) ■ Vad^a - (z x VaA" 1 ^) • V a 5Z a 
= [(z X VaC e ab ) ■ VaA” 1 - (z X VAa"'^) ■ V„] SZ a . 

Similarly for SA b 
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with k± = k ± n/2. We also have that 


d 2 k 


= e 


(2tt)2 
in-(x a +x 6 )/2 


U,Q, = /^ Mn x k) f 4 - i ) C"(k) e'< k +">*.- k x. 


l 


A: 2 


n z 


/ 


^^(nxlc_)(-l-i)c«( k _) e ,k - 


,ik-(x a -x b ) 


(E.6b) 


and similarly, 


CS, = / ^ S2 ■ (k x n) (i - i) C*(k) e ik --‘< k +”) k * 


— _ e in -( x a+Xb)/2 


/ 


d 2 k 

(2^F 


1 1 


z • (n x k+) ( p- - ] C e (k+) e 1 ' 


,ik-(x a -x b ) 


77^ 


(E.6c) 


From (E.lb) we then obtain that 


C (h) (k) = 


z-( nxk -) C ,e (k_)-z-(nxk+) ) £ e ( k +) 




cr 


+ 2 + ^2+ i ( w k + — w k_ — w n) 


,in-(x a +x 6 )/2 


F(k_) - F( k+) 


0 in-(x a +x 6 )/2 


(E.7) 


with a = s + ito ! n . Further, 


K(5C n ) = 


= -V • 

= -V • 

= -V • 


^(VaA-' + VtA^SC 

Z Jx a =X(, 

i X (v a A-> + V^ 1 ) / [F (k _) - F(k +) ] 


gin-(x a +xfr)/2gik-(x a —Xfr) 


x a =x 6 


d 2 k z /ik_ 
(2^)22 X 



F(k_)-F(k+)] e in x , 


(E.8) 


and since under the transformation k —k we have that k + — k_, k\ k 2 _ 

and F( k+) —F( k_), (E.8) becomes 

^ = " V 7(W 2 X (^ - ^) F ^ X ■ (E.9) 
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Using V • (z x pe in ' x ) = — iz • (n x p) e in ' x , we then obtain, 


K(6C n ) = 


/ 


d 2 k . , , / 1 1 ' 

w" (nxl) 


F(k_)e i] 


- S2 I 


d 2 k 


|n X k|2 [kj fc 2 ) (fc 2 n 2 


eQ(k) 


(2vr) 2 a + 2 + v 2h (k 2h + k 2h ) + i (w k+n - - w n ) 2(1 + u 2h k 2h ) 


(E.10) 


with k s = k+n and k s = |k s |. The last equality was obtained with the substitution 
k k — n/2. For h = 0 (E.10) gives (3.9). Finally, using (3.3a) we obtain the 
dispersion relation for the stability of the homogeneous equilibrium: 

2h= r d 2 k |nxk|2 (fcf ~p) (fc 2 ~^) _ Q(k) 

1 2,1,1 £ J (27I-) 2 a + 2 + V 2h {k 2h + k 2h ) + i (w k +n - - w„) 2(1 + u 2h k 2h ) 

(E.ll) 


Equation (E.ll) can be written in terms of the real and imaginary part of the 
eigenvalue s as: 


s r = -(1 + V 2h r? h ) + e Re[/(cr)] , (E.12a) 

Si = ~uJ n + £ Im[/(er)] . (E.12b) 


The real part of the eddy feedback contributes to the growth rate of the mean 
flow and the imaginary part determines the departure of the phase speed of the 
mean flow from the Rossby wave frequency. For /? ^> 1 the marginally unstable 
eigenfunctions have cq = —cj n as it can be readily shown that fi is at most of 
0(1) and therefore produces only a small correction to the Rossby phase speed 
which is of 0(f3). 

Equation (E.ll) is solved numerically for a for a given n, /3, e and Q( k). 
However, in some special cases a can be solved in closed form. Such an example 
is when /3 = 0 and with z/2 h — 0. Then (E.ll) takes the form: 


(s + l)(s + 


2)=ej 


d 2 k |k x n| 2 ( k 2 — k 2 )(k 2 — n 2 ) Q(k) 


(2vr) 


k^k 2 n 2 


= eJ 


(E.13) 


and the eigenvalues are s — 


—3/2 ± (1/4 + eJ) 1/2 . 


The value of the integral, J, 
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depends only on the wavevector n and the forcing spectrum. For the ring forcing 
spectrum (3.11) by using the definitions of Fig. 3.4 we have 


k x n| 2 = n 2 cos 2 6 , k 2 — 1 + n 2 + 2n sin 6 , 


(E.14) 


and therefore J takes the form: 



o 


= ^n 2 (l — n 2 ) cos(2 ip) , for n < 1 . 


(E.15) 


For isotropic forcing J — 0 and there is no instability since cr — — 1 or d = -2. 
When we have instability when £ > 32/[/icos(2(/?)] > 0, which occurs either 

when /I > 0 and 0 < (^ < 45° or when p < 0 and 45° < (f < 90°. 

E.2 Derivation of expression for f r for the ring forcing spectrum 

In this section we will derive the expression (3.13) for f r = Re [/(0)] for the ring 
forcing spectrum (3.11). For V2h — 0, cr — 0 and for spectrum (3.11) we have that 


/<«=/ 


kdkdO |k x n| 2 (k 2 — k 2 )(k 2 — n 2 ) 
(2t) 2 k A k^n 2 2 + i (c^k+n — u n — c^k) 


Q(0 — ip) 8{k — 1) 



(E.16) 


Using the definitions of Fig. 3.4 we have: 


cj n = —/3siiup/n , — — /3cos(6 — cp) , c^k+n = —0 


n sin cp + cos(0 — ip) 


1 + n 2 + 2n sin 0 

(E.17) 


and together with (E.14) we get 



o 



(E.18) 


154 














with F(0, n) defined by (3.14) and 

Vq(9, n) = 2(1 + n 2 + 2nsin0) , (E.19a) 

V 2 (9 , n) = (1 + n 2 + 2nsin0) sin p/n + n 2 cos {6 — p) + nsin (2 9 — p) , 

(E.19b) 

NiO, n) = — n ^1 — n 2 ^j cos 2 0(sin0 + n/2) Q{9 — p) . (E.19c) 


Note that /3V 2 is 


(3V 2 (9,n) = k 2 (u> k+n -6c; n -cj k ) , (E.20) 

and therefore vanishing of V 2 occurs when the resonant condition is satisfied: 

— ^k+n • (E.21) 

Also note that the F defined by (3.14) remains unchanged when the angle 
p is shifted by 180° (p -A 180° + p) or when there is a simultaneous shift of 
p -A 180° — p and 6 -A 180° — 6. As a result, it suffices to only consider cases 
with 0 < p < 90°. 


E.3 Asymptotic expression for the eigenvalue s r at high supercriti¬ 
cality 


At high supercriticality, i.e., as e -A 00 , the maximal growth rate, s r , of the 
large-scale structure with wavevector n scales with y/e while the frequency of this 
eigenstructure, Si, asymptotes to values Si ~ — oj n . 

Specifically, s r asymptotes to: 


= £ I 


d 2 k 

WT 2 


|nxk| u^ 


M ( 1 M <?( k ) 

k 2 ) \k 2 n 2 ) 2 


(E.22) 


This asymptotic expression for the growth rate and phase speed of the large- 
scale structure is useful for tracing the maximal growth rates as a function of 
supercriticality using Newton’s iterations. 

The asymptotic growth rates depend only on the forcing distribution and for 
the forcing spectrum (3.11) are shown in Fig. E.l for fi > 0 and (i < 0. It can 
be shown that the asymptotic growth rate vanishes for exactly isotropic forcing. 
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Asymptotically the growth rates do not 
flow, r m (cf. chapter 5). 


(a) M = i 



Tlx 


depend on the damping rate of the mean 


(b) n = -l 



Figure E.l: The £ +oo asymptotic maximal growth rate s r scaled by as a function 
of the wavenumbers (n x ,n y ) of the S3T eigenfunction. Shown are contours for s r > 0 and the 
zero contour is marked with a thick solid line. The asymptotic growth rate is independent of (3 
and dissipation and depends only the forcing spectrum. Shown are the asymptotic growth rates 
for forcing (3.11) for (a) fi — 1 and (b) \i = — 1. For /i > 0 maximal instability occurs for jet 
structures, while for /x < 0 maximal instability occurs for non-zonal structures. 
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F 


Asymptotic expressions for eddy feedback 


In this Appendix we calculate in closed form asymptotic expressions for the eddy 
feedback induced by a mean flow perturbation with wavevector n, in the cases 
/? <C 1 and (3 1. 

F.l Case f3 < 1 

When /3 1 and for n satisfying f3/n <C 1, we expand J r (0,n) = F{0, n) + 

F(180° + 0, n) in (3.13) in powers of f3. Since T is a function of /3 2 we have the 
expansion: 


T = To + f3 2 + 0((3 4 ) , 

with J ~2 = \ dppF ■ The leading order term is: 


(F.l) 



(F.2) 


due to the property C/(180° + 6) — G{0). Positive values of Fq indicate that 
the stochastically forced waves with phase lines inclined at angle 6 with respect 
to the wavevector n (cf. Fig. 3.4) induce upgradient vorticity fluxes to a mean 
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flow with wavenumber n when /3 = 0. Given that n < 1 and Q > 0, is 
positive for any forcing distribution, only in the sector shown in Fig. 3.5a in 
which 4 sin 2 6 < 1 + n 2 . Specifically, in the absence of /? all waves with \6\ < 30° 
reinforce mean flows with n < 1. Note that the condition 4sin 2 0 < 1 + n 2 is 
also the necessary condition for modulational instability of a Rossby wave with 
wavevector components (cos #, sin 6) to any mean flow (zonal or non-zonal) of 
total wavenumber n for /? <C 1 (Gill, 1974). 


The total vorticity flux feedback f r for Q(6 — p) = 1 + fi cos [2(6 — p)\ is at 
leading order: 

fr = (l - n 2 ) cos(2 ip) + 0(/3 2 ) , (F.3) 

which is proportional to the anisotropy factor, p. The maximum feedback factor 
is in this case 

•,max — ’ 

and is achieved for mean flows with n — l/\/2. This maximum is achieved for 
zonal jets (p = 0°) if p > 0 and for meridional jets (p = 90°) if p < 0. This 
implies that for f3 1 the first structures to become unstable are zonal jets if 
/i > 0 and meridional jets if fi < 0, as shown in Fig. 3.12c. 

For isotropic forcing (fi — 0), the leading order term is zero and f r depends 
quadratically on /?: 



fr = P 2 — 2 + cos(2 p) + 0(f3 4 ) for n < 1 , 
64 L -l 


(F.5) 


producing upgradient fluxes for n < 1. Note that for the delta function ring 
forcing / Q 27r J^d^ is discontinuous at n — 1, with positive values for n — 1~ and 
negative values for n — 1 + . The accuracy of these asymptotic expressions is 
shown in Fig. F.l. The maximum feedback factor, shown in Fig. 3.12a, is 


L 


r, max 


3f3 2 
~64 ’ 


(F.6) 


and is attained by zonal jets (<p = 0°) with wavenumber n —> 1~ as (3 —¥ 0, a 
result that was previously derived by Srinivasan and Young (2012). The accuracy 
of (F.4) and (F.6) extends to /3 « 0.1, as shown in Fig. 3.12a. 
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Figure F.l: Feedback factor f r for a non-zonal perturbation with n — 0.4751 and ip = 10° 
(which belongs in region A of Fig. 3.8a) (solid lines) in the case of a forcing covariance with (a) 
li — 0 and (b) p — 1/4. Also shown are asymptotic expressions for /3 1 ((F.6) in (a) and (F.4) 

in (b)) and the resonant contribution (3.17) for p 1 (dash-dot). For ft 1, expression (F.12) 
is also plotted (dashed). It can be seen that only (3.17) can captures the /3 _1//2 decrease of f r . 


F.2 Case /3 > 1 


When f3 1, we write (E.18) in the form: 

27r 

fr = jp, with I = j F x (d, n) d d , (F.7) 

0 

where 

(F - 8) 

and x = 1//5- When V 2 ~ 0(1) for all angles 0, then the feedback factor is 
f r ~ 0(/3~ 2 ). However, if V 2 ~ 0(f3~ l ) for some angle 0, then as we will show in 
this Appendix, f r decays as 0(f3 ~ 1 ) or as (9(/3 -1 / 2 ). This is illustrated in Fig. F.l 
showing the feedback factor f r as a function of /3 in cases in which V 2 vanishes. 

V 2 can have at most 4 roots, 0° < 0j < 360° (j — 1,2, 3,4), for any given 
(n, <£>). At these angles the resonance condition (E.21) is satisfied. To calculate 
asymptotic approximations to the integral /, we split the range of integration to 
a small range close to the roots of V 2 for which we have resonance, I^ R \ and to a 
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range away from the roots of V 2 , /( NR ); 


N r 

" Oj-SO 
r 

0j+50 

r 

II 

M 

J F x {e, n)de+ 

J F x (9 : n)d9 

3 = 1 

_8j- 1+<56» 

Oj-se 


j(NR) j(R) 


(F.9) 


where 7V r is the total number of the roots of V 2 and 9q = 9j\f r . Asymptotic 
approximations to the integral over the two ranges are then found separately 
using a proper rescaling for the regions close to the roots of V 2 (cf. Hinch (1991)). 

When the distance between two consecutive roots is \9j — 9j-i\ > y/x, as in the 
examples shown in Figs. 3.8c,e, then the dominant contribution to the integral 
comes from the 0(x) regions close to the roots Qj, since F x (9, n) close to 9j is 
approximately a Lorentzian of half-width 0(x)- Therefore, choosing the range 89 
close to the roots to be yfx <C 89 <C 1, Taylor expanding F x (9, n) close to 9j and 
rescaling 9 — 9j + x u we obtain: 


r( R ) _ 

3 


so/x 



—50/x 


Mj Vqj d u 

v h + ^ 2 


+ ^(x“ 3 ), 


(F.10) 


where V' 2 = BqV 2 and the subscript j denotes the value at 9j. In the limit 
89/x 00 we obtain: 

(R) _ 1 7T Mj ( , 

~ xMj ’ (F11) 

and as a result, the resonant contribution produces the asymptotic approximation: 


f(R) = 1 V 77 

fr i' 


(F.12) 


However, special attention should be given to the case in which two consecutive 
roots are close to each other. When 1 9j - Oj- il ~ O(Jx) then V' 2J ~ 0(y/x) 
and fr R ^ scales as instead of 1//3 for f3 1. Indeed, when F x is double 

peaked, as in Fig. 3.8d, the dominant contribution comes from the whole range 
between the two resonant angles which are a distance 0(y/x) apart. The proper 
scaling for the angles close to 9j is therefore 9 — 9j + yfx u - Taylor expanding the 
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denominator under this scaling we obtain: 


x 2 t> 1 +vI = x 2 ^oj+x^u 2 +x 3 / 2 v' 2 j v'l jU 3 +x 2 (-v'Q + \v2M3) ^+o{x b ' 2 ), 


(F.13) 

where V 2 = cIqqDi and T> 2 = c^XY When V 2 j ~ O(^x) all the terms in 
(F.13) are 0(% 2 ) an( i writing V 2 j = y%d( n > 6j) = ^Jx dj . where d is of 0(1), the 
leading order resonant contribution is: 



se/Vx 

3/2 f _ A/} Vqj du _ 

J T>1 ■ + df u 2 + djPjU 3 + Ipju 4 

-se/Vx ,J J 4 3 


+ 0{ X ~ l ) , 


(F.14) 


where pj = I n th e limit 86/yjx oo the integral can be evaluated 

from the residues from two of the four poles of the integrand. The two poles 
are at u = —dj/pj ± l^jl 1 / 2 sgn (pj) e ±lw ^ 2 , where \zj\ = Vqj |pj| _1 (/v- + 4) 1 / 2 , 
wj = arctan (2 / kj) and Kj = d?-V^ x -\pj\~ x is an increasing function of the distance 
between the two roots of T> 2 . Therefore: 


;T = * 


-3/2 _ 


nMjVj 

Y/W 1/2 


+ o(x 


- 1 \ 


(F.15) 


and 


f( r) _ J_ V'' 1 r( R ) _ V" 7 P^/l r li 


(F. 16 ) 


which is exactly (3.17). The factor 1/2 in (F.16) arises because the range of 
integration includes both angles and (F.15) must be divided by 2 , in order to 
avoid double counting. The resonant response is proportional to 


77 = 2 (k 2 + 4) 3//4 esc 


"1 

2 


arctan( 2 /ft) , 


(F.17) 


which is always positive, because k > 0 as Do > 0. The factor 77 is shown as a 
function of k, (which is a rough measure of the distance between the roots) in 
Fig. F.2. We observe that the maximum value is attained at hi = 2/\/3 ~ 1.16, 
that is when the roots are at a distance 0(x 1 ^ 2 ) apart. Note also that by taking 
the limit of the resonant angles being away from each other, that is by taking 
the limit k 1, 77 ~ 2 / yfd and (F.16) reduces to (F. 12 ). Consequently, (F.16) 
is a valid asymptotic expression regardless of the distance between the roots 6j. 
The accuracy of (F.12) and (F.16) in comparison with the numerically obtained 
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integral is shown in Fig. F.l. 

The sign of the resonant contribution depends only on the sign of J\f. From (E.19c) 
we see that J\f > 0 when sin 6 > —n/2 for n < 1; this region is highlighted with 
light shading in Fig. F.3. It should be noted that for the important case of 
zonal jet perturbations (cp = 0°) the resonant contribution is exactly zero because 
J\fj — 0, as shown in Fig. F.3a. The asymptotic behavior of the feedback factor 
for this case is found from the non-resonant part of the integral. Expanding in 
this case the integrand for y 1, we obtain to leading order: 

fr » / r (NR) = (1 - n 2 )( 2 + ii)r 2 + <3(/r 4 ) , (F.18) 


with the maximum feedback gain 


fr ,max — (2 + /i)/3 2 + 0(/3 4 ) , (F.19) 


occurring for n —> 0. 

Consider now non-zonal perturbations (cp / 0°). There is a large region in the 
(n,(p) plane (region D in Fig. 3.8a) in which V 2 has no roots and f r = 0(/3~ 2 ). 
For larger values of n (region B in Fig. 3.8a), and for any given ip, V 2 = 0 for 
exactly two Qj that satisfy the inequality sin Qj < —n/2. Consequently, Afj < 0 
and the resonant contribution from these roots is negative. For even larger values 
of n (regions A and C in Fig. 3.8a), V 2 has exactly 4 roots. Only two of the 
roots in region A produce positive resonant contributions. Note also that region 
A extends to ip < 60° and (p > 120 0 . 1 

The maximum response, which is 0(/3 -1 / 2 ), arises in region A close to the 
curve separating regions A and C where n ~ 1.16. While the roots of V 2 are 
independent of /?, the location and the size of the region of maximum response 
depends on (3 through the dependence of k on /?. However, as (3 increases this 
dependence is weak and as (3 -A 00 the maximum response occurs in a narrow 

x It can be shown that fluxes from the resonant contributions for n < 1 are necessarily 
downgradient (negative) for 60° < ip < 120°. Proof: A positive contribution is produced when 
the T >2 — 0 curve enters into the J\f > 0, highlighted with light grey in Fig. F.3. There are 4 
roots of T >2 on the unit circle n — 1 (on which also M = 0), at angles: 0 — 210°, 270°, 330° and 
6 — 90° + 2 Lp (marked with A, B , C and D respectively). The V 2 = 0 curve can cross the curve 
AOC , which separates positive from negative AT , only at points A and C, since V 2 — 0 only at 
these points on AOC. Therefore, the T >2 — 0 curve can enter the AT > 0 region i) through D , if 
it lies outside the arc ABC , and/or ii) through A, C. However, for 60° < ip < 120° point D lies 
within the arc ABC and moreover, the gradient 2 at points A and C is oriented in such way 
that does not allow the V 2 — 0 curve to enter Af > 0, as d n V 2 < 0 and deV 2 < 0 (doT >2 > 0) at 
point A (point C). 
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Figure F. 2 : The factor 77 = 2 (k 2 + 4 ) -3 / 4 esc [| arctan ( 2 /«)] as a function of ft that is a 
measure of the distance between two consecutive resonant angles. The maximum value of 77 
marked with an open circle (and consequently of the feedback gain that is proportional to 77 ) is 
77 = 3 3//4 /2 « 1.14 and it is achieved at ft = 2 /v/3 ~ 1.16. Also shown is the asymptote 77 = 2 /y^ 
that 77 follows for ft 1 (dashed). This suggests that the resonant contribution is maximum 
when the two roots are very close to each other (ft ~ 1 ) but not on top of each other (k, <C 1 ). 


region near n ~ 0.5 and ~ 10°, marked with a star in Fig. 3.8a. The width of 
this region decreases with /3, making it exceedingly hard to locate for large /3, 
and the asymptotic approach of (n,(p) to (0.5,10°) is shown in Figs. 3.12b,c. 
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Figure F. 3 : Zero contours of £> 2 ( 0 , n) for (a) zonal jet perturbations (p = 0°), (b) non-zonal 
perturbations with p = 15° and (c) non-zonal perturbations with p = 75° in a (0, n) polar plot. 
Shaded areas mark n < 1. Light shade corresponds to (0,n) satisfying sin0 > —n/2 for which 
we have positive resonant contributions ( J\f > 0), while dark areas correspond to sin0 < —n/2 
for which we have negative resonant contributions ( J\[ < 0). Points of intersection of the T >2 — 0 
curve with the unit circle are marked with A, B , C, D. The radial grid interval is An = 0.25. 
The curve — 0 does not enter the N > 0 area for 60° < p < 120°. 
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G 

Formal equivalence between the S3T 
instability of a homogeneous equilibrium with 

the modulational instability of a 
corresponding basic flow 


In this Appendix we demonstrate the formal equivalence between the modulational 
instability (MI) of any solution of the barotropic equation, which may be in general 
time dependent but has stationary power spectrum, with the S3T instability of the 
homogeneous state with the same power spectrum. Consider a solution ^(x, £), 
with vorticity Qq — A^, of the inviscid and unforced nonlinear barotropic 
equation (2.1) with time-independent power spectrum. Because J(^g,Cg) — 0, 
(g satisfies the equation 


dt( G = £ (h) Cc , 


(G.l) 
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with = z • {(3 x V) A x . Linear perturbations S( to this solution evolve 
according to the equation: 


d t s( = £S(, 


(G. 2 ) 


where 

C — -ug • V + (Aug) • VA " 1 + z • (/3 x V) A ” 1 = C' G + £ (h) , (G.3) 

C' G £00 

is the time-dependent linear operator about (g that has been decomposed into a 
spatially homogeneous operator, that governs the evolution of £q and the 
inhomogeneous operator CJ G that depends on (q. The hydrodynamic instability 
of £g is ascertained when the largest Lyapunov exponent of (G. 2 ) is positive. 

We proceed with the study of the MI by decomposing the perturbation into 
a mean 5Z = (5Q and deviations from the mean S( f = — 5Z , where ( • ) 

is an averaging operation. The averaging operation in MI is projection to the 
eigenstructure with wavenumber n, which is orthogonal to £< 3 , because only 
orthogonal eigenstructures to (q could become unstable. With this averaging 
operator (£< 3 } = 0, and therefore C G = Cg> whereas the perturbations has a non¬ 
zero mean, 5Z, and a deviation and is expressed as 5( = 6Z + 5£'. For example, if 
g is a sum of Rossby waves as in (4.9) the perturbation field from Bloch’s theorem 
comprises Fourier components with wavenumbers n, nip j, n±2pj, n±3pj,... 
for all the p j. In this case 8Z is a plane wave with wavenumber n and S(' 
comprises the remaining Fourier components. With these definitions (G.2) is 
equivalently written as: 

d t (SZ + S(') = C’qSZ + £ (h) <5C' + C G S(' + £ (h) SZ , (G.4) 

where Cl G is primed in order to stress that the operator linearly depends on the 
deviation quantity ( G . Equation (G.4) is then separated to form an equivalent 
system of equations for the evolution of the mean perturbation, 5Z, and the 
deviation perturbation, 5(': 

d t 8Z = & h Hz + (£ g 8£) , (G.5a) 

dt 8(' = £ <ll >8(' + a G 8Z + a G 8C! - <£'g ^C'> • (G.5b) 
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The stability equation (G.2) and the stability equations (G.5) for 8Z and 8C/ are 
equivalent. 

In MI studies usually the term CJ G 8C/ — (CJ G 8(') in (G.5b) is neglected and the 
stability of the following simpler system is studied: 

d t 5Z = 6 h) 5Z + (C g 5C) , (G.6a) 

d t SC' = £ {h) 5(' + CqSZ . (G.6b) 

For example, if is in the form of (4.9) the neglected term comprises waves 
with wavevectors n ± 2p j5 n ± 3p J 5 ... and the truncated system (G.6) allows 
only interaction between the primary finite amplitude waves p J? the perturbation 
n and the waves n =b pj. If ( G is a single wave p (as in MI studies), (G.6) is 
referred to as the 4 mode truncation or “4MT” system since it comprises only 
modes p, n and nip. 

However, instead of studying the MI stability of 8Z and 8(' using the approxi¬ 
mate equations (G.6), we can equivalently study the stability of 8Z and 

5C{x a , X fe , t) = ( C o( x a, t) S('(x b , t ) + Cg( x &> t) S C '(x a , t) } 

= ( C G,a + Cg, 6 <Ka ) • (G-7) 

With these definitions we obtain from (G.l) and (G.6b) the evolution equation 
for 8C: 

d t 8C = ^(9«CG,a) S(b + (dt(G,b) S(a + C G,a ^tKb) + C G,b (^Ca)) 

= ( 4^ + 4^ ) SC + ^ Cg,o 4?,6 $ Z b + Cg,6 4 G,a S Z a ) • (G.8) 

We note from the definition of C G (cf. (G.3)) that: 

C G 8Z = — (zx V^ G ) ■ V8Z + (z x VCg) • V<Stf 
= (z x V8Z ) • - (z x • VCg 

= (A 8XJ) • 0VV4) - (S\J) • (V&) = SACg, (G.9) 

where 5U = z x Vd'h is the velocity field associated with SZ and 8A = —8 U-V + 
[(A <7U)-V] A -1 is the operator that also appears in (3.2b). As a result, (G.8) 
becomes: 

d t sc = ( 4 h ) + 4 h) ) SC + (, 5A a + 5At) C G , (G. 10 ) 
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where C G = {C'g a^Gb)- Returning now to (G.6a) we note that (C'qSQ') = 1Z(5C), 
where 7 Z(SC) is defined in (2.8), as: 


7 Z(8C) = -V • 
= -V • 


| x ( v.a - 1 + VfcA," 1 ) (CG,a K + C kb <) 

|zx((V^)5C' + (V# / )Cg)} 


= -V-(u^^' + ^u / C^) 

= ( - • V 5C! + (Au^ • V<ty') = {C' G 6?) . (G.ll) 


Consequently, the MI of CJ G in the approximation (G.6) is equivalently deter¬ 
mined from the stability of the system: 


dt sz = C {h) 5Z + 1Z(5C) , (G.12a) 

dtSC = ( 4 h) + 4 h) ) 5C + (dA a + SA b ) C G , (G.12b) 


which is identical to equations (3.2) that determine the S3T stability of the 
homogeneous equilibrium with zero mean flow, U e = 0, and equilibrium covariance 
C e = C G under the ergodic assumption that ensemble averages are equal to 
averages under operation ( • ). 

For example, consider the nonlinear solution 

27T 

^(x, t) = J a{6 ) cos(p • x — (jj p t) d 6 , (G.13) 

o 

with wavevectors p = (cos 0, sin 9) on the unit circle (p = 1) and take f3 = (0,/3). 
Expanding the plane waves into cylindrical waves: 


gi [(x+/3t) cos 6+y sin 0] 


+oo 

m =—oo 


(G.14) 


with R 2 — ( x+/3t) 2 +y 2 , </> = arctan y/(x + /3t ) 
of the first kind, this can be shown to be the non-dispersive structure 


and Jm the m-th Bessel function 


^(x + f3t,y) = Re 


+oo 


^ ^ Ifm Jm{R) 


Jimcf) 


(G.15) 


propagating westward with velocity /?, where = / Q 27r a{6) e irn6 d6. The re 


D —imO . 
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suits in this Appendix show that the MI of the propagating structure (G.13) 
in the approximation (G.6) is equivalent to the S3T instability of the homoge¬ 
neous equilibrium with covariance C e prescribed by power spectrum (7 e (k) = 
(27t) 2 \a(6)\ 2 5{k — 1). Note that this S3T equilibrium is also an exact homo¬ 
geneous statistical equilibrium of the nonlinear barotropic equations without 
approximation. 
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Specification of the forcing structures used in 

chapter 5 


In chapter 5 three spatial structure of stochastic forcing are used in the inves¬ 
tigation of the correspondence among S3T, QL and NL dynamics: a forcing 
with narrow isotropic ring spectrum (IRFn), a forcing with wide isotropic ring 
spectrum (IRFw) and a forcing with non-isotropic spectrum (NIF). 

For the IRFn we take the forcing spectrum 


Qk 


c if \k — kf \ < 5kf 
0 if | k — kf | > 5kf or k x = 0 . 


(H.l) 


The constant c is chosen so that Qk satisfies 


E Q k _ 

2k 2 ” 


k x fki 


(H.2) 


and therefore, according to (A.20), the energy injection rate by Qk is a unit. For 
the narrow ring forcing we choose kf — 14 and 5kf = 1. 
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For the IRFw we take the forcing spectrum to be 


Qk 


j cexp 

-(k-k f ) 2 /(2 5k})_ 

\o 



if k x ± 0 
if k x = 0 


(H.3) 


Again c is chosen so that (H.2) is satisfied. For the wide ring forcing we choose 
kf — 14 and 5kf = 8 . 

For NIF we force the zonal wavenumbers k x — 1,..., with power: 

Qk = c kx d 2 e~ k l d2 , (H.4) 

with constants c^ x chosen in such manner so that for a fixed zonal wavenumber 
k x -, 

^2 (kZ + k*) = W k ’ (H - 5) 

so that all zonal wavenumbers k x inject the same energy input rate and the total 
injection rate by Qk is unity. For the anisotropic forcing we force wavenumbers 
k — 1,..., 14 and d = 1/5. 

Comparison of the forcing spectra structure (H.l), (H.3) and (H.4) as well all 
forcing realizations that they induce are shown in Fig. 5.1. 
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I 


Determination of inhomogeneous zonal jet 
S3T equilibrium solutions using Newton’s 


iteration 


In this appendix we present a method for determining stationary equilibrium 
solutions (u e (y),C e (x a — £ 6 ,y a >y&)) of the S3T system (2.13). Remember that 
for solutions of the form (u(y, t), C(x a — y a , 7/5, t)^j and for the usual ori¬ 
entation of (3 = (0,/?) the S3T system (2.13) collapses to the simpler S3Tz 
system (2.21) and therefore it is sufficient to determine stationary equilibrium 
solutions of the S3Tz system, or equivalently (cf. Appendix C.2.1), stationary 
solutions (j7 e , Cf,..., C e Nk ) of (C.22) that satisfy the matrix equations: 



We determine the equilibrium solution satisfying (I.la) through Newton’s 
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iterations. However instead of iterating both the U and C^ x towards the solution, 
which is computationally very costly, for each flow iteration U we determine 
the Ck x that satisfy (I.lb) and treat with this understanding the Ck x as linear 
functions of J7, writing Ck x = Ck x (U). This approach has the disadvantage that 
we constrain the iteration to the manifold of hydro dynamically stable flows for 
which there exist physical steady covariances that satisfy (I.lb). Under these 
conditions we can obtain the equilibrium state by solving for U e the linear 
equation: 


N k 

G (U e ) = 2Re 
— 1 


vecd (ik x A^ 1 C kx (U e )) 


- U e = 0 . 


( 1 . 2 ) 


In order to solve (1.2) with Newton’s method we start the iteration by selecting a 
hydrodynamically stable flow C7°. If G (u^^j — 0 no iteration is needed. If not, 
assume that the equilibrium U is nearby, and hence to first order it must satisfy: 


0 = Gi (17) « Gi (l7(°)) + 


'dG 

dU 




(1.3) 


where dG/dU is is the Jacobian matrix of G with elements [dG/dU] - — dGi/dUj , 
Ui and Gi are the N y elements of U and G respectively and subscript denotes 
that the matrix elements are evaluated at U — . The iteration is continued 

by taking as the new iterant the U that satisfies (1.3), i.e., 


[/(D 


jy(°)_ 



g(uW) , 

U ( o) 


(1.4) 


where ( dG/dU ) 1 is inverse of the Jacobian matrix dG/dU. We approximate 
the elements of dG/dU\ u=U (o) with 


dGi 

Gi | 

(uW + hej) -Gi | 

[u^-hej) 

dUj 

u( 0 ) 

2 h 



(1.5) 


where ej is the unit-vector in the j -th direction with elements [ej]i = Sij and h > 0 
is sufficiently small. After U ^ is calculated from (1.4) we repeat the iteration. If 
the initial guess U° is close to the equilibrium U e the iteration converges rapidly 
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to the equilibrium solution. We consider the iterations converged if 


Si 

G J 

(u^y 

)i 2 





< 1(T 14 . 


( 1 . 6 ) 
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J 

Stability of inhomogeneous S3T equilibrium 

solutions 


Consider a solution (Z. C ) of the S3T system: 


dtZ + J (®, Z + (3 ■ x) = K(C) - Z , 
9tC a b = A (U) + A(U)] C a b + e Qab 


(J.la) 

(J.lb) 


with A defined in (2.5). Linear perturbations (5Z,SC) about this solution obey: 


dt 5Z = U(U) 5Z + n{5C) , 
d t SC ab = [A(U) + A b { U)] 5C ab + (<5A + 5A b ) C ab , 


(J.2a) 

(J.2b) 


with 5A — w4(U + SXJ) — w4(U). The stability of (Z,C) is determined by the 
top Lyapunov exponent of (J.2), which is calculated numerically using the power 
method. We initiate the integration of (J.2) with a normalized perturbation 
state (SZ, SC), using with any norm. At every time-step, /i, we calculate the 


175 




perturbation state growth, 


log (\\(SZ(jh),5C(jh))\\\ 

A (jh) = -A--- L , j = 1 , 2 ,..., (J.3) 

and then renormalize the perturbation state before moving to the next time-step. 
After sufficient time-steps the growth A converges to the top Lyapunov exponent 
and the state ( SZ , SC) to the first Lyapunov vector. 

The time-integration of the discretized (J.2) proceeds as described in Ap¬ 
pendix C.2.2 and system (J.2) takes the form: 

5Z = A(U) SZ + R(<fC) , (J.4a) 

SC = A(U) SC + SC [A(U)] t + <5A C + C (<5A) T , (J.4b) 

where R is defined in (C.33), A is defined in (C.30) and S A = A(U + SJJ) — A(U). 
In (J.4) SZ is an (A^A^)-column vectors while JC is an (N x N y ) x (N x N y ) matrix. 

J.l Stability of zonal jet equilibria 

The homogeneity in x of zonal jet equilibria enables us to seek S3T eigenfunctions 
in the form 

SZ(x,t) = e 1HxX J dn y a Ux (n y , t) e myV , (J.5a) 

5C(x a ,x b ,t) = ^ Xa+Xb )/2 j dkx SC nx (k x ,y a ,y b ,t)e ik ^- x ^ . (J.5b) 

i.e., the perturbation mean flow is a single harmonic in x with wavenumber 
n x and the inhomogeneous part of the perturbation covariance is also a single 
harmonic with the same wavenumber. We will show that by using (J.5) the S3T 
perturbation system (J.2) results in a smaller system of decoupled equations for 
Fourier components <T nx (n y ,£) and 5C nx (k x ,y a ,y b ,t). 

In an infinite domain the wavenumbers n y and k x assume continuous 
values. However, system (J.2b) is solved in a channel of dimension L x x L y 
with periodic boundary conditions. For a box with L x = L y = 27t, the periodic 
boundary conditions on SZ impose that n x and n y take only integer values and 
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the periodicity on 5C imposes that 


k 


X 


integer , for n x even , 

< 


half-integer, for n x odd . 


(J.6) 


(A number m is called half-integer when m + 1/2 G Z, which implies that 
m = (2 k + l)/2 for k G Z.) Therefore, after redefining k x k x + n x j 2 so that k x 
assumes integer values for any value of n x , the eigenfunction (J.5) takes the form: 


SZ(x,t) = e in * x J2 a n*,n y (t)e in yy , (J.7a) 

n y 

5C(x a ,x fe ,i) = , (J.7b) 

k X 

where notation (C.4) is used in the summation. Notice that we have chosen 
to define SCk x ,n x as the Fourier coefficient that corresponds to e 1 ( k x-n x /2)(x a -x b ) ^ 
The equilibrium covariance C e is also expanded as in Appendix C.2.1, 


c e (x a ,x 6 ) = Y J C e kx {ya,y b y kx{Xa ~ Xb) ■ (J.8a) 

k x 

We want to write now the S3T perturbation system (J.2) in terms of the 
Fourier coefficients OL nx ,n y (t) and 5Ck x ,n x (Va , £). Consider first a single harmonic 
of (J.7a) with vorticity: SZ n = —n 2 e mx . Then: 


S^ n ^SJJ n ,V5Z n = — l/n 2 , — (z x V)/n 2 , in 6Z n , (J.9a) 


and 


at Bo «C = c ta -^)/ 2 5: , (J.lOa) 

k x 

d Xb 8C = £ -i (fc x - n,) 8C kxtnx (y a , y b ) , (J.lOb) 

k x 

A a 5C = ^^(4 - *£) 5C kx>nx (y a ,y b )j<- k *- n ’W x °- x ») , (J.lOc) 

k x 

A b 8C = ^(x a +x b )/2 [d 2 yb - (k x - n x ) 2 ] 8C kx , nx (y a ,y b ) , 


(J.lOd) 
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imply that: 


-(U e a d Xa + U§d Xb )SC nx = 

= [-i k x U e a SC kx , nx (y a ,y b ) 


+i(fc x n x )U b 8C kxtUx (y a , y b ) 


J(k x -n x /2)(x a -x b ) 


(J.lla) 


~(SU a d Xa + 6U b d Xh )C e = 

= ~ E [ 8U - ik * Cl (?/«> Vb) + SU b (-i k x ) C e kx (y a , y b ) 


Jk x (Xa—Xb) 


— in y e 


in x (x a +x b )/2 


Y {[i ik * 


8E nyA ) C e kx (y aj y b ) 


+ 


[C e kx (y a , Vb) SE ny>b j 


J(k x +n x /2)(x a -x b ) 


J(k x —n x /2)(x a -x b )'\ 


Jn x (x a +X b )/2 


Y { 


n x) 8E n a C k —n (ija- lib) 


+ 


Ck x (y a ,yb ) (-i k x SE„ VJ b^ 


a \{k x -n x /2){x a -x b ) 


(J.llb) 


with the notation SE Uy = e 1Uyy . We see that each component of (J.2b) has as a 
common factor the term e m x(x a +x b )/2. p or an y mean fj ow perturbation (J.7a) we 
can rewrite (J.2b) as a system of matrix equations for the perturbation covariances 
matrices 6C kx ,n x , with elements [5Ck x ,n x {t)] ab = SC kx ,n x (ya,yb,t), and equilibrium 
covariance matrices C kx with elements [C kx ] ab — C^ x (y ai yb)- Consequently, for 
the usual orientation (3 = (0,/3), each 5C k x ,n x satisfies: 


~ ^k x ^k x ,n x + 3Clk x ,n x (\-fcJ 


+ ^ a n x ,n y ^k x -n x ,n x ,n y ^k x -n x + ^k x (^^-k x ,n x ,n y ^j 


with: 



(J.12) 

A% x = —ik x U e - i k x (, 61 - U^) A^ 1 - 1 , 

(J.13a) 

8&k x ,n x ,n y = ~8E ny (- \-n y k x \ + in x Dy ) fl + n A^ ^ , 

(J.13b) 
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and 5E Uy = diag(e m ^). 

Because SC and C e are symmetric to the exchange x a X 5 and further because 
C e is real, we have that 


c i, = = (Q,)' , (J.14a) 

^-k x ,n x ~ (^Ckx+nxi'nx') • (J.14b) 

From (J.14a) we see that is only necessary to determine C% x for non-negative k x . 
In general, if C e has non-zero C| only for k x = 1 ,..., then for a fixed n x we 
only need to solve for 2 (7V& + n x ) + 1 perturbation covariance matrices SCk x , nsc , 
since perturbation covariances with \k x \ > N^ + n x are necessarily zero because 
they do not couple with any of the non-zero C% x . In practice however, we only 
solve for the first + 2 n x matrices SCk x ,n x and deduce from them the rest N^ + 1 
using the symmetry (J.14b). Note that for n x — 0 perturbations, in the case 
that C% x=0 / Owe see that we need to solve for + 1 perturbation covariances 
SCjz X: 0 , and not for N\ as it is claimed above. However, it can be proven that 
SCk x =0,n x =0 — 0, even if C| 0 ^ 0 and therefore only N perturbation covariances 
are non-zero. 

Turning now to the mean flow perturbation equation (J.2a) we want to express 
the Reynolds stress divergence associated with perturbation covariance SC in 
terms of the covariances SCk x ,n x - At point x m the Reynolds stress divergence is 


K(5C) 


{ dy a A^+^A ,- 1 


J Vb L 


,-(d Xa ^~ l + d Xb ^^))sc ab 


J x a =x b 


(J.15) 


Using (J.5) and expressing everything in term of the matrices SCk Xl n x we get 

d x [\(d ya A- 1 + d yi A^)sc ah 


v 2 

= d x {~ 


Xa=Xb 


{ ■I TViu — \~ Tlx 

l E (^a .- 1 + «wV) 

k x =-(N k +n {c ) 


Jn x (x a +x b )l 2 sc kx , nx (y a , y b ) eKfe-nx/ 2 )^-^)' 


Xa=X(, 
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1 , 


N k +n x 


— Q x ) _ e i n x (x a +x b )/2 ^2 


D V^k x ^k Xi n x 


— -m T e 


1 Ti X Xm 


k x =-(N k +n x ) 

+dCk x ,n x (^1 

N k +n x 

E 


J(k x -n x /2)(x a -x b ) 


ab 


x a =x fe - 


fea; = —(iVfc+n*) 


fi^k x ,n x + 3Ck x ,n x (^k x -n x ) 


V) 


(J.16) 


Similarly, 


-dyU^A-'+d^A^SCab _ \ 

t x a— x 6 ) 


= —^ vecd 


A^-b^cc 


XI ^ k x ,n x + k x ,n x k x -n x ) I 

_k x = (A^^+rz-x) 

(J.17) 


From (J.16) and (J.17) it can be seen that the e 1UxX dependance factors out and 
therefore the Reynolds stress divergence takes the form 7 Z(SC) = e 1TlxX TZ Ux (SC nx ). 
After discretization 7 Z nx (5C nx ) is approximated by the column vector SR Ux : 

N k +n x 

SR n. = 2 £ 

k x =-(N k +n x ) 

+ D^vecd 


(For a homogeneous S3T equilibrium, i.e., in the absence of mean flow U e = 0, the 
Reynolds stress Rn x (SC nx ) becomes proportional to e myV or SR Ux is proportional 
to vecd[JE n J.) 

Note that the operator D^A^ 1 appearing in (J.16) is ill-defined for k x — 0, 
since then A^ 1 is the non-invertible operator (D” 1 ) 2 . In this case we calculate 
D^A^ 1 as the pseudoinverse of D y using its SVD decomposition and then by 
removing the zero singular values of D y . 


i n x vecd 


^y^k x ^k x ,n x + ^k x ,n x (Py^k x -n x ^) 


^ k x ,n x + k x ,n x (y^kx-rix) [ n x)\ 


(J.18) 
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(J.19a) 


After all these consideration the discrete version of (J.2) is: 


±5Z nx = tf n JZ nx +6R nx , 

d 'p 

~^^k x ,n x ~ ^k x ^k x ,n x + kx,n x (A n x —k x ) 


Ny/ 2-1 

+ a n x ,n y 

n y =—N y /2 


3^k x -n x ,n x ,n y ^k x —n x 


+ Cf 


(^A —k x ,n x ,riy^ 


for h x — 1 •> • • • •> Nk + ^ti x , 

(J.19b) 


where 5Z Ux is the A^-column vector with elements [SZ Ux (t)] a = Y Uy a n x ,n y (t) e myVa 
This system has a state of 2N y + 2(Nk+2n x )Ny real variables. For N x = N y = 128, 
Nk = 15 and n x — 2 this gives an order of 400-fold decrease in the dimension of 
the state variable compared to the full S3T perturbation system (J.4). 

An example demonstrating the convergence of the growth A to the top Lyapunov 
exponent using (J.19) is shown in Fig. 6.lie. 
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